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CHAPTER  1 


INTRODUCTION 


Military  aircraft  of  the  future  are  envisaged  to  be  multirole  vehicles.  This  dictates 
the  requirement  for  multimode  integration  of  airframe  and  propulsion  systems  over  a 
wide  flight  envelope.  In  most  present  designs,  the  multimode  requirement  is  being  met 
by  varying  aircraft  configurations  or  relying  on  active  controls  technology  to 
accommodate  response  changes  due  to  varying  flight  conditions.  It  is  becoming  evident 
that  in  order  to  meet  this  variable  capability  much  of  it  must  be  met  by  the  propulsion 
system. 

Advanced  aircraft  turbine  engines  have  been  the  focus  of  extensive  research  in 
recent  years  under  the  Integrated  High  Performance  Turbine  Engine  Technology 
(IHPTET)  program.  The  IHPTET  program  is  a  national  initiative  sponsored  by  the 
Department  of  Defense  (DoD)  and  the  National  Aeronautics  and  Space  Administration 
(NASA)  seeking  to  improve  the  overall  performance,  reliability  and  affordability  of 
advanced  gas  turbine  engines  (Brown,  1990).  This  manifests  itself  in  the  need  for  smaller 
engines  that  produce  the  same  or  higher  levels  of  thrust  while  having  better  fuel  economy 
and  costing  less  to  maintain. 
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IHPTET  is  a  multiphase  program  with  very  specific  goals  to  be  accomplished 
during  each  phase.  In  order  to  achieve  these  goals  the  turbine  engine  was  broken  down 
into  its  major  component  areas  and  objectives  were  established  for  each  component  that 
would  then  be  related  to  or  influence  the  overall  goals  of  the  program  (Petty  and 
Henderson,  1987).  For  the  propulsion  control  system  these  objectives  are  weight 
reduction  and  design  margin  reduction.  The  weight  reduction  goal  has  been  addressed  by 
material-for-material  substitution  in  the  major  control  components  such  as  fuel  pumps 
and  actuators.  The  design  margin  goal  looks  more  at  the  overall  impact  that 
improvements  in  the  control  system  architecture  and  logic  can  have  on  the  engine  system. 
This  second  goal  is  the  one  being  addressed  by  the  work  described  herein. 

1.1  Review  of  Relevant  Research 

The  application  of  multivariable  time-  and  frequency-domain  control  design 
techniques  to  aircraft  turbine  engine  control  problems  has  been  widely  published  in  the 
literature  [see  (Merrill,  et  al.,  1984)  for  an  extensive  list  of  references].  Many  of  the 
published  solutions,  however,  suffer  from  excessive  controller  order  and  complexity, 
necessitating  the  use  of  complex  control  logic  to  implement  the  control  law.  The 
development  of  multivariable  frequency-domain  control  design  techniques  such  as 
Rosenbrock’s  (1969)  Inverse  Nyquist  Array  (IN A)  promised  the  use  of  relatively 
inexpensive  analog  hardware  or  simple  digital  control  software  to  implement  the  control 
law.  This  has  not  proven  to  be  the  case  in  that  many  of  the  published  solutions  are  only 
valid  at  a  single  operating  point  in  the  flight  envelope,  necessitating  the  use  of  digital 
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hardware  and  gain  schedules  to  implement  a  different  control  law  for  each  point  in  the 
flight  envelope  (Merrill,  et  al.,  1984;  Polley,  et  al.,  1988).  This  has  the  undesirable  side 
effect  of  requiring  extensive  test  cell  and  installed  engine-controller  testing,  which  is  very 
expensive.  The  applicability  of  frequency-domain  control  design  techniques  to  this  type 
of  nonlinear  control  problem,  however,  has  been  demonstrated  (Leininger,  1981). 

The  development  of  Quantitative  Feedback  Theory  (QFT)  to  handle  uncertain 
multivariable  control  problems  has  enjoyed  tremendous  success  pver  the  past  decade, 
particularly  in  the  area  of  controller  design  for  plants  that  may  be  accurately  modeled  by 
connected  sets  of  parametrically  uncertain  differential  equations  (Horowitz,  1991). 
Indeed,  QFT  based  techniques  have  been  applied  to  turbofan  engine  control  problems  in 
the  past,  such  as  the  work  of  Yau  et  al.  (1994),  however,  an  appeal  to  plant  high  gain 
stabilizability  was  resorted  to.  Due  to  the  inherent  time  delay  behavior  of 
turbomechanical  systems,  this  is  not  a  very  realistic  assumption.  Recently,  the  Inverse 
Nyquist  Array  design  technique  (Rosenbrock,  1970)  was  modified  by  Nwokah  and 
Nordgren  (1993)  using  techniques  from  QFT  and  Hco  control  theory  (Zames,  1981)  to 
handle  multivariable  control  problems  for  uncertain  plant  sets  with  both  parametric  and 
unstructured  uncertainty.  This  decentralized  design  technique,  termed  the  Quantitative 

Nyquist  Array  (QNA),  offers  considerable  advantage  over  design  techniques  that  treat 

) 

model  uncertainty  in  a  single  manner. 

The  technique  of  p-synthesis  is  rooted  in  Hoo  control  theory  (Zames,  1981)  and 
has  been  pursued  extensively  by  Owen  and  Zames  (1981);  Doyle,  et  al.,  (1992);  and 
Chaing  and  Safanov,  (1988).  p-synthesis  seeks  to  eliminate  some  of  the  performance 
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weaknesses  of  the  standard  Hoo  theory  such  that  desired  performance  measures  are 
achieved  for  all  plant  set  elements.  This  theory  has  been  expanded  to  handle  plants  with 
both  parametric  and  unstructured  uncertainty  models.  Quantitative  Feedback  Theory  and 
p-synthesis  offer,  perhaps,  the  greatest  promise  in  the  development  of  a  single 
multivariable  robust  control  law  that  is  capable  of  delivering  desired  engine  performance 
levels  throughout  the  flight  envelope. 

Multimodal  engine  control  was  first  looked  at  under  the  Air  Force  funded 
Intelligent  Engine  Control  (IEC)  (Adibhatla,  et  al.,  1992)  program  and  was  subsequently 
followed  by  studies  conducted  by  NASA-Lewis  Research  Center.  From  these  studies  it  is 
obvious  that  conventional  control  techniques  would  not  satisfy  all  the  needs  for  the 
future.  Therefore,  a  new  paradigm  was  needed  to  meet  future  requirements.  This  is  when 
the  technique  of  model-based  control  was  introduced  by  Adibhatla,  et  al.  (1992)  and  Qi 
and  Maccallum  (1993). 

1.2  Significance  of  the  Work 

The  main  objective  of  turbine  engine  control  systems  is  to  adjust  the  system  for 
rapidly  required  thrust  changes,  while  maintaining  stable  operation  of  the  compressor  and 
keeping  the  turbine  inlet  temperature  within  safety  limits.  Usually  specific  engine 
performance  rating  points  are  generally  defined  as  steady-state  design  goals.  In  addition 
to  regulating  the  engine  for  steady-state  operation  the  control  must  constrain  the  operation 
within  various  aerodynamic,  thermodynamic,  and  mechanical  or  structural  design  limits. 
Typical  aerodynamic  limits  include  fan  and  compressor  stability  and  maximum  allowable 
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engine  inlet  airflow  distortion.  Thermodynamic  limits  include  minimum  combustor  or 
augmentor  fuel/air  ratios  for  stable  combustion,  minimum  fuel  flow  for  ignition  and 
maximum  fuel/air  ratios.  Mechanical  or  structural  limits  consist  of  maximum  allowable 
rotor  speeds,  maximum  allowable  nozzle  flap  temperatures,  maximum  combustor  case 
pressure,  rotor  creep  limits,  and  maximum  average  and  peak  values  of  turbine  blade  metal 
temperature.  Steady-state  performance  rating  points  at  sea  level  static  and  altitude  along 
with  transient  response  characteristics  are  specified  in  “General  Specification  for  Turbojet 
and  Turbofan  Engines”  (MIL-E-5007D). 

The  first  turbine  engines  used  hydromechanical  control  systems  to  meet  all  the 
requirements  levied  on  the  propulsion  system.  However,  the  need  for  increased 
functionality,  reliability  and  maintainability  led  to  the  introduction  of  digital  engine 
controls  in  the  early  1980’s.  According  to  Skira  and  Agnello  (1991)  today’s  turbine 
engine  control  systems  are  greatly  improved  over  the  first  electronic  control  systems. 
However,  these  first  digital  electronic  controls  did  little  more  than  emulate  the 
hydromechanical  control  they  replaced.  Since  that  time  much  work  has  been  done  on 
improving  these  systems  by  adding  redundancy  and  evolving  the  components  that  go  into 
them  with  little  change  in  the  control  strategy.  The  primary  emphasis  has  been  on  weight 
reduction  with  very  little  attention  given  to  the  power  provided  by  today’s  advanced 
computers. 

Protection  of  the  engine  hardware  is  one  of  the  primary  responsibilities  of  the 
control  system.  This  consists  of  avoiding  all  the  limits  listed  above.  In  order  to  avoid 
reaching  these  limits,  margins  (safety  factors)  are  established  that  consist  of  a  worst  case 
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stackup  of  undesirable  effects  that  could  cause  the  engine  to  exceed  one  or  more  of  the 
limits.  In  the  case  of  stall  margin,  this  consists  of  effects  due  to  but  is  not  limited  to 
values  associated  with:  inlet  distortion,  engine  transients,  Reynolds  number  effects, 
engine-to-engine  variations  and  control  tolerances.  Historically,  the  amount  set  aside  for 
these  design  margins  has  not  changed. 

The  control  system  can  not  create  performance  that  the  machinery  does  not 
already  have  built  into  it.  However,  performance  is  designed  out  of  a  system  in  order  to 
build  in  large  safety  factors.  Historically  these  margins  are  constant  even  though  there 
has  been  significant  technological  advancements.  This  margin  stackup  is  a  series  of 
worst  case  events  that  could  occur  for  each  effect.  It  rarely  occurs  that  the  maximum 
effect  in  each  parameter  occurs  at  the  same  time.  Therefore,  if  the  margin  stackup  is  used 
in  a  smart  way  we  can  regain  some  of  the  lost  performance  while  protecting  the  engine 
and  not  sacrificing  the  safety  or  integrity  of  the  machinery. 

In  order  to  illustrate  these  concepts  consider  Figure  1.1.  In  this  plot,  thrust  ratio  is 
plotted  as  function  of  temperature  ratio.  Temperature  ratio  is  defined  as  the  turbine  inlet 
temperature  limit  divided  by  the  turbine  inlet  temperature  at  stoichiometric  temperatures. 
This  axis  then  reflects  in  general  the  turbine  engine  material  limitations  that  exist  over 
time.  Thrust  ratio  is  then  defined  as  the  thrust  produced  at  a  given  stall  margin  and 
temperature  ratio  divided  by  the  thrust  at  no  stall  margin  and  stoichiometric  temperatures. 
The  line  labeled  stall  line  is  the  line  of  zero  margin  and  the  one  labeled  operating  line  is  a 
line  of  constant  design  stall  margin.  The  area  between  these  two  curves  is  the 
performance  that  is  lost  due  to  constant  design  stall  margins.  The  ellipses  represent 
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engine  generations  or  level  of  engine  technology  as  time  progresses.  The  wedge  labeled 
recaptured  performance  is  shown  this  way  to  illustrate  that  understanding  of  the  effects  of 
the  parameters  that  make  up  the  margin  stackup  are  not  fully  understood  now.  As  we 
gain  further  understanding  the  amount  of  performance  that  can  be  regained  will  increase. 


Figure  1.1  Recaptured  Performance  Potential 

In  order  to  recapture  this  loss  performance  a  better  understanding  of  the  processes 
that  influences  each  factor  in  the  margin  stackup  is  necessary  and  a  means  of  determining 
the  conditions  where  each  factor  is  dominant  needs  to  be  found.  Once  this  has  been 
accomplished  then  the  control  can  offer  additional  performance  by  using  up  the  margin  or 
providing  engine  protection  when  the  operating  conditions  require  it. 
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1.3  Conventional  Versus  Model-Based  Control 

During  any  engine  development  process  a  great  deal  of  time  is  spent  creating  a 
detailed  nonlinear  simulation  of  the  engine.  The  control  engineer  usually  starts  his  work 
after  the  simulation  is  created.  First  the  nonlinear  model  is  linearized.  Then  the  linear 
model  responses  are  checked  against  the  nonlinear  model  responses.  Finally,  the  control 
design  proceeds  with  the  linear  models.  It  is  at  this  point  in  the  design  process  that 
significant  improvement  can  be  made  if  the  problem  is  rethought  slightly. 

This  nonlinear  model  that  is  developed  contains  all  the  knowledge  and  expertise 
available  for  that  engine  at  that  time.  But,  it  is  never  used  for  anything  more  than 
evaluation  of  proposed  changes  to  control  logic  or  engine  hardware.  The  question  is,  “Is 
there  a  way  that  the  knowledge  contained  within  the  nonlinear  simulation  can  be  used  to 
improve  the  functionality  of  the  engine?”  The  answer  is  yes,  through  the  control 
methodology  termed  model-based  control.  Now  consider  the  differences  between 
conventional  and  model-based  control. 

The  conventional  control  system  is  shown  in  Figure  1.2.  In  this  type  of 
architecture,  sensed  parameters  such  as  fan  speed  and  pressure  ratio  are  used  as  feedback 
signals  to  generate  error  signals  that  the  control  uses  to  drive  actuators  to  their  correct 
position.  These  sensed  parameters  are  used  because  they  infer  what  the  parameters  of 
interest  (thrust  and  stall  margin)  should  be.  Because  this  inference  is  never  exact,  large 
margins  (safety  factors)  are  established  in  order  to  protect  the  engine. 


8 


Reference 


Figure  1.2  Conventional  Control  System 


The  model-based  control  architecture  is  shown  in  Figure  1.3.  The  key  to  this 
architecture  is  that  a  model  (a  computer  program  representation  of  the  plant)  is  included 
in  the  feedback  path.  This  model  can  compute  variables  of  interest,  which  can  not  be 
measured  or  are  not  easily  measured.  These  computed  values  are  then  used  to  generate 
error  signals  that  the  controller  uses  to  drive  the  actuators  to  their  desired  positions.  The 
main  advantage  of  having  a  resident  model  is  that  the  control  can  be  done  on  variables  of 
interest  (i.e.,  thrust  and  stall  margin).  Then  the  margins  can  be  managed  more  efficiently 
thereby  allowing  for  an  increase  in  performance. 

The  key  to  model-based  control  is  the  tracking  filter.  Its  sole  purpose  is  to 
guarantee  that  the  model  is  matching  the  plant.  This  becomes  somewhat  of  a  problem 
when  there  are  only  a  limited  number  of  sensed  parameters  available  to  match.  In 
addition,  the  parameters  to  update  can  become  enormous.  In  general,  the  tracking  filter  is 
square,  the  same  number  of  inputs  and  outputs.  However,  current  research  is  directed 
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toward  the  development  of  methods  for  generating  nonsquare  tracking  filters.  A  working 
nonsquare  tracking  filter  has  not  yet  been  found.  The  tracking  filter  is  expected  to  work 
well  in  whatever  mode  the  engine  is  operating  in.  This  also  presents  a  problem  in  that  the 
updates  tend  to  distort  the  model  in  such  a  way  as  to  optimize  the  parameter  of  interest. 
This  parameter  does  not  have  to  be  the  same  for  every  mode  of  operation.  Then  a 
tracking  filter  that  works  well  for  every  mode  of  operation,  termed  a  multipurpose 
tracking  filter,  is  needed.  However,  to  date  no  such  filter  exists.  There  are  solutions  to 
this  problem  such  as  multiple  copies  of  the  model  each  with  their  own  tracking  filter  for 
each  mode  of  operation.  But  these  solutions  are  not  very  viable  for  aircraft  engine  control 
applications  due  to  the  lack  of  physical  space  to  put  additional  memory  in  for  additional 
models  and  tracking  filters. 


Figure  1.3  Model-Based  Control  System 
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1.4  Other  Control  Methodologies  that  Use  Internal  Models 

There  are  two  main  control  design  techniques  that  exist  in  today’s  literature  that 
make  use  of  a  model  in  the  control:  Model-Reference  Adaptive  Control  (MRAC)  and 
Internal-Model  Control  (IMC).  There  are  key  differences  between  Model-Based  Control 
(MBC)  described  in  this  work  and  these  existing  techniques.  However,  one  of  the  main 
differences  in  every  case  is  how  the  error  signal,  generated  by  comparing  the  outputs  of 
the  plant  and  model,  is  used. 

Model-Reference  Adaptive  Control  (MRAC)  is  one  of  the  main  approaches  to 
adaptive  control.  The  MRAC  system  shown  in  Figure  1.4  consists  of  a  plant  (P),  a 
reference  model  (M),  a  feedback  regulator  (controller)  (C)  and  an  adaptation  system.  The 
reference  model  contains  the  desired  performance  and  will  express  the  desired  output  to 
the  external  command.  The  feedback  regulator  has  adjustable  parameters  that  are  varied 
through  the  adaptation  system  based  on  the  error  between  the  reference  model  and  the 
plant.  The  plant  outputs  are  also  used  for  traditional  feedback  to  compare  with  the 
reference  input  to  generate  a  signal  for  the  feedback  regulator.  The  MRAC  approach  was 
developed  to  handle  the  case  where  the  specifications  are  given  in  terms  of  a  desired 
response  model  (Astrom  and  Wittenmark,  1989). 

Internal-Model  Control  (IMC)  was  initially  applied  to  solve  problems  in  the 
process  control  industry,  and  is  related  to  disturbance  accomodation  in  the  system,  where 
empirical  or  dynamical  models  of  the  processes  under  control  are  available  (Morari  and 
Zafiriou,  1989).  The  IMC  system  shown  in  Figure  1.5  consists  of  a  plant  (P),  a  plant 
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model  (PM)  and  a  feedback  regulator  (C).  In  this  approach,  the  error  signal  generated  by 
comparing  the  plant  and  plant  model  outputs  is  used  as  the  feedback  signal. 


Figure  1.4  Model-Reference  Adaptive  Control  Structure 


Figure  1.5  Internal-Model  Control  Structure 


The  Model-Based  Control  (MBC)  approach  described  in  Section  1.3,  which  is 
different  from  MRAC  and  IMC,  was  developed  to  solve  the  problem  of  controlling  plant 
outputs  that  axe  not  directly  measurable.  The  key  to  this  architecture  is  that  the  plant 
model  variables  are  changing  due  to  the  adaptation  error  signal  (the  difference  between 
the  plant  an  the  plant  model  outputs,  when  subjected  to  the  same  input)  fed  into  the 
adaptation  system  (the  tracking  filter  described  in.  Chapter  3).  When  this  error  signal  is 
below  a  certain  threshold  the  plant  model  computes  the  desired  feedback  variables  that 
are  compared  to  the  standard  reference  inputs,  to  generate  the  reference  error  signal  which 
is  fed  to  the  feedback  regulator  for  tracking  response. 

1.5  Overview 

In  this  work  a  concise  and  workable  methodology  for  the  design  of  a  robust 
multivariable  model-based  control  system  will  be  presented.  The  key  areas  to  be 
addressed  are  the  design  of  the  tracking  filter,  and  the  design  of  a  robust  control  law  using 
unmeasurable  feedback  variables.  While  model-based  control  techniques  have  been 
discussed,  they  have  never  been  pursued  in  the  context  of  robustness. 

The  remainder  of  this  document  is  arranged  in  the  following  manner.  A  complete 
description  of  the  nonlinear  engine  model  developed  for  this  effort  and  its  linearization  is 
given  in  Chapter  2.  The  tracking  filter  is  considered  in  Chapter  3.  In  Chapter  4  the 
robust  multivariable  control  design  problem  is  discussed  along  with  H-matrix  theory.  A 
detailed  design  problem  is  presented  in  Chapter  5.  Chapter  6  summarizes  the 
contributions  of  this  research,  and  discusses  some  future  research  topics.  Appendix  A 
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contains  a  well-known  methodology  for  finding  the  transition  matrix  required  to  convert 
an  observable  system  into  the  observable  canonical  form.  Appendix  B  contains  proofs  of 
important  theorems  given  throughout  this  dissertation.  The  Computer  routines  used  to 
perform  the  work  presented  herein  are  given  in  Appendix  C.  Finally,  in  Appendix  D  the 
normalized  linear  models  are  given. 
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CHAPTER  2 


TURBINE  ENGINE  FUNDAMENTALS  AND  MODELING 

2.1  Introduction 

The  propulsion  system  of  choice  for  a  wide  variety  of  land,  sea,  and  air  vehicles  is 
the  turbine  engine.  Because  of  varied  operational  and  performance  requirements,  future 
aircraft  systems  will  impose  unprecedented  challenges  upon  the  system  designer. 
Advanced  military  aircraft  contemplated  for  the  future  and  beyond  will  likely  be  designed 
as  multimission  aircraft  with  both  air-to-air  and  air-to-ground  capabilities.  These  same 
systems  will  possess  many  of  the  following  operational  requirements:  Short  Take-Off  and 
Vertical  Landing  (STOVL)  capabilities,  automatic  threat  evasion,  terrain  avoidance/threat 
avoidance,  and  automatic  weapon  delivery.  In  addition,  reliability  and  maintainability 
requirements  will  be  much  higher  for  future  systems  (Adibhatla,  et  al.,  1992).  Therefore, 
an  understanding  of  the  basic  principles  of  turbine  engine  design  and  modeling  is 
paramount  to  understand  the  work  presented  herein. 

2.2  Engine  Fundamentals 

The  primary  purpose  of  the  turbine  engine  (shown  in  Figure  2.1)  is  to  impart  a 
change  in  momentum  to  a  mass  of  fluid.  This  change  in  momentum  is  equivalent  to  an 
acceleration  of  a  working  fluid  producing  an  external  force  (Thrust)  on  the  system.  In  a 
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turbojet  this  is  accomplished  by  bringing  air  into  the  inlet,  compressing  the  mass  of  air, 
mixing  fuel  with  the  high  pressure  air  in  the  combustor  section  there  igniting  and  burning 
it,  the  fluid  is  expanded  in  the  turbine  section  driving  the  turbines  which  in  turn  power  the 
compressors.  The  fluid  is  further  expanded  through  the  nozzle  section  to  a  high  velocity 
(conversion  of  pressure  and  thermal  energy  into  kinetic  energy)  thus  increasing  the 
momentum  of  the  fluid  and  producing  thrust. 


Figure  2.1  J85  Nonafterburning  Engine 


In  addition  to  acting  as  the  propellant  fluid,  the  air  acts  as  the  working  fluid  in  a 
thermodynamic  process.  This  process  can  be  represented  thermodynamically  by  the  Air 
Standard  Brayton  Cycle  as  shown  in  Figure  2.2  for  the  turbojet.  The  processes  that  make 
up  this  cycle  are  given  in  Table  2. 1 . 
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Figure  2.2  The  Brayton  Cycle 


Table  2.1  Brayton  Cycle  Processes 


Leg 

process  description 

1-2 

Reversible,  adiabatic  (isentropic)  compression  between  minimum  and 
maximum  pressures 

2-3 

Heat  addition  at  constant  maximum  pressure 

Reversible,  adiabatic  (isentropic)  expansion  between  maximum  and  minimum 
pressures 

4-1 

Heat  rejection  at  constant  minimum  pressure 
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The  behavior  of  turbine  components  is  normally  represented  on  "maps".  The 
performance  maps  for  compressor  and  turbines  normally  consist  of  the  following 
parameters:  component  pressure  ratio,  corrected  mass  flow  rate,  a  corrected  speed 
parameter,  and  adiabatic  efficiency.  A  typical  compressor  map  is  shown  in  Figure  2.3. 
The  curve  marked  stall  line  represents  a  limit  on  compressor  performance.  Operation 
below  this  line  is  essential  for  satisfactory  engine  performance.  Steady  operation  above 
this  line  is  impossible  and  entering  this  region  even  momentarily  is  hazardous  to  the 
engine  structural  integrity. 


Figure  2.3  Compressor  Map 
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An  example  turbine  map  is  shown  in  Figure  2.4.  The  main  difference  is  that  on 
this  map  it  is  the  reciprocal  pressure  ratio  plotted  as  a  function  of  corrected  mass  flow  rate 
and  corrected  speed.  Since  the  turbine  is  usually  choked,  other  methods  of  modeling  the 
turbine  can  be  used  just  as  effectively  if  a  map  is  not  available.  One  such  method  is  using 
the  continuity  equation  and  considering  the  turbine  pressure  distribution  as  a  series  of 
pressure  drops  across  the  elements.  This  is  the  method  used  for  this  work. 


Figure  2.4  Turbine  Map 

The  engine  community  has  developed  a  systematic  notation  to  facilitate  the 
manipulation  of  equations  and  to  describe  at  what  point  along  the  gas  path  that  is  being 
discussed.  The  work  herein  will  follow  the  standardized  gas  turbine  engine  station 
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identification  and  nomenclature  system  recommended  by  SAE,  Inc.,  Aerospace 
Recommended  Practice  (ARP  75  5 A).  Each  station  location  refers  to  a  particular  engine 
location.  These  locations  are  identified  in  Table  2.2.  See  Figure  2.5  for  a  schematic 
representation. 


Table  2.2  Station  Designation  and  Associated  Locations 


Station 

Location 

0 

Far  Upstream 

1 

Inlet  Entry 

2 

Compressor  Entry 

3 

Compressor  Exit 

4 

Combustor  Discharge 

5 

Last  Turbine  Discharge 

6 

Available  for  Mixer,  Afterburner,  etc. 

7 

Engine/Exhaust  Nozzle  Interface 

8 

Exhaust  Nozzle  Throat 

9 

Exhaust  Nozzle  Discharge 
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2.3  Engine  Modeling 

A  simulation  or  mathematical  representation  of  gas  turbine  engines  is  fairly 
common  and  has  been  reported  by  several  authors:  Szuch  and  Seldner  (1975),  French 
(1982),  Mattingly,  et.  al.  (1987),  Oates  (1988),  Seldner,  et.  al.  (1971),  Sellers  and  Terren 
(1974)  and  Szuch  (1974).  For  this  work  a  simulation  of  a  single-spool  turbojet  engine 
was  developed  using  thermodynamic,  aerodynamic,  mechanical,  or  empirical 
relationships  of  each  of  the  major  components.  The  model  was  adjusted  to  match  actual 
test  data  of  a  turboj  et.  This  matching  was  accomplished  through  efficiency  adjustments. 

Maps  were  used  to  establish  the  airflow  and  thermodynamic  performance  of  the 
compressor.  The  map  was  changed  into  lookup  tables  of  temperature  ratio  and  corrected 
airflow  for  conditions  of  corrected  speed  and  pressure  ratio.  Performance  of  other 
components  was  established  to  agree  with  the  compressor  data  and  engine  test  data.  Due 
to  the  proprietary  nature  of  the  maps  they  can  not  be  printed  in  this  document. 
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2.3.1  Compressor  Efficiency.  The  compressor  efficiency  at  partial  power  was 
modified  to  a  lower  value  in  order  to  produce  a  fuel  flow  that  is  more  representative  of 
the  value  obtained  during  actual  engine  testing.  The  efficiency  was  modified  so  that,  for 
the  idle  speed  used  in  the  simulation  the  efficiency  would  be  75  percent  as  opposed  to  the 
84  percent  obtained  from  the  map.  The  efficiency  matched  well  at  military  speed. 
Therefore,  the  efficiency  was  blended  from  idle  (68%  corrected  speed)  to  military  speed 
(100%  corrected  speed)  where  no  correction  was  needed.  The  efficiency  correction  is 
accommodated  in  the  simulation  by  dividing  the  temperature  rise,  T3  -  T2,  computed  from 
the  map  by  the  correction  factor. 

2.3.2  Compressor  Torque.  The  compressor  torque  was  obtained  by  equating  the  heat 
energy  per  unit  time  of  the  compressor  gas  to  the  work  per  unit  time  required  to  drive  the 
compressor.  The  work  per  unit  time  defined  by 


ft.-lb. 

minute 


(r,)(N)2* 


[2.1] 


The  power  required  based  on  the  heat  consideration  is  defined  as 


h  =  cpclhc(T3-T2)  [2.2] 

Equating  equations  [2.1]  and  [2.2]  in  the  proper  units  and  solving  for  the  torque,  results  in 
the  following  relationship. 


KmccPc(T3-T2) 
T‘= - N - 


[2.3] 
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In  addition,  equation  [2.3]  must  be  modified  by  the  efficiency  correction  mentioned 
above  by  dividing  the  temperature  rise  by  the  correction. 

2.3.3  Combustor  Equations.  The  combustor  performance  is  defined  by  the  balance  of 
the  energy  added  by  the  fuel  equated  to  the  energy  required  to  produce  the  temperature 
rise  in  the  combustor.  The  energy  produced  by  the  fuel  is  considered  to  be, 

hf=mf(Q)qb  [2.4] 

The  energy  required  to  produce  the  (T4  -  T3)  temperature  rise  is  defined  as, 

K  =  cpbI^14(Xi  ~T3)  [2.5] 

Equating  equations  [2.4]  and  [2.5]  then  solving  for  temperature  rise  results  in  the 
following  formula, 

(T,-T3)=S!M-  [2.6] 

cP,m4 

In  order  to  increase  the  steady-state  idle  fuel  flow  in  the  simulation  to  make  it  agree  with 
the  engine  test  data,  the  efficiency  of  the  combustor  is  varied  with  speed  similar  to  the 
method  used  for  compressor  efficiency. 


23 


2.3.4  Turbine  Gas  Flow.  Ideally  a  turbine  map  would  have  represented  the  turbine; 
however,  none  was  available.  The  flow  equation  was  established  from  known 
information.  The  continuity  equation  was  the  basic  equation  used.  The  turbine  pressure 
distribution  is  considered  to  be  the  result  of  a  series  of  pressure  drops  that  represents  the 
two  rotor  and  two  nozzle  diaphragms.  A  split  of  45  percent,  55  percent  was  considered 
for  the  two  stages.  The  pressure  distribution  across  each  stage  was  considered  to  be  60 
percent  for  the  diaphragm  to  40  percent  for  the  rotor.  Using  this  relationship,  the  drop 
across  the  assembly  was  represented  as: 


1st  diaphragm 

(ps/p4)(.^,=(Pj/p4).. 

1st  rotor 

(ps/p)r,'-,=(p!/p<)» 

2nd  diaphragm 

(vp.r'^ =(p,/p.)' 

2nd  rotor 

5! 

'o 

U* 

"o 

II 

5 

o 

The  following  equation  was  derived  from  the  continuity  of  flow: 


'  out 


[2.7] 


where  k  is  the  ratio  of  specific  heats,  C  is  a  proportionality  constant,  and  A  is  an  area. 
This  equation  expresses  the  flow  through  each  of  the  turbine  restrictions.  The  flow 
through  the  first  diaphragm  is  used  to  determine  the  mass  flow  through  the  turbine  with 
Pout/Pin  being  equal  to  (P5/P4)027. 
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2.3.5  Turbine  Torque.  The  turbine  torque  equation  also  had  to  be  determined  without 
the  aid  of  turbine  maps.  The  equation  was  established  from  the  impulse  and  momentum 
equation  for  the  rotor.  The  Torque  on  the  turbine  was  expressed  as: 

Torque  =  f(m5 ,  Vgas ,  Vrotor )  [2.8] 

The  velocity  of  the  gas  was  assumed  to  be  approximately  twice  the  velocity  of  the  rotor. 
The  velocity  of  the  gas  was  determined  at  each  point  through  the  turbine  in  order  to 
determine  an  average  velocity.  This  average  velocity  and  the  following  equation  were 
then  used  to  calculate  gas  velocity. 


gas 


K, 

k 

l- 

fPsT 

■  0.073 

i 

j . 

[2.9] 


These  assumptions  were  then  combined  to  determine  turbine  torque. 

2.3.6  Turbine  Temperature  Drop.  The  equation  for  the  turbine  temperature  drop  was 
determined  in  much  the  same  manner  as  the  compressor  equation.  The  coefficient  of 
specific  heat,  cPi,  was  considered  to  be  0.26  at  the  temperature  in  the  turbine.  The 
equation  describing  the  turbine  temperature  drop  is  given  by: 

(T,-T,)  =  ^  [2.10] 

Cp,m4 
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2.3.7  Torque  Balances.  The  speed  of  the  engine  is  determined  by  integrating  the 
acceleration.  The  acceleration  was  established  from  the  relationship: 

(t.-t) 

N  =  — - cJ-  [2.11] 

KJ 

The  rotor  inertia  was  estimated  from  experimental  test  data. 

2.3.8  Nozzle.  The  flow  through  the  exhaust  nozzle  was  determined  from  the  continuity 
equation  for  the  nonchoked  condition  similar  to  equation  [2.7].  A  limit  of  0.534  was 
placed  on  the  (P0/P7)  value  used  in  this  equation  to  account  for  the  choked  condition. 

2.3.9  Airflow  Balances.  In  order  to  relate  the  airflow  through  the  compressor,  turbine 
and  nozzle,  airflow  balance  equations  were  required.  The  transient  effect  was  taken  into 
account  by  considering  the  compressibility.  The  general  equation  used  for  generating  the 
pressures  using  compressibility  is. 


YA 

rtc 


[2.12] 


This  equation  was  used  for  airflow  balances  within  the  combustor  and  the  tailpipe.  In 
addition,  the  turbine  exit  pressure  (P5)  was  determined  from  the  nozzle  inlet  pressure  (P7). 
It  is  assumed  that  there  is  a  two-percent  loss  in  the  tail  section.  In  a  similar  maimer 
compressor  exit  pressure  (P3)  is  determined  from  turbine  inlet  pressure  (P4).  It  is  assumed 


that  there  is  a  four-  percent  loss  in  the  combustor. 


All  these  equations  were  then  used  to  build  a  nonlinear  engine  model  within  the 
Simulink  environment  (Mathworks,  Inc.,  1994).  Figure  2.6  show  the  block  diagram 
representation  of  the  engine.  In  this  model  each  of  the  major  components  are  separately 
modeled  the  communication  loop  and  the  input  and  output  management  blocks  handle  the 
passing  of  variables  between  the  components.  This  model  will  be  used  in  subsequent 
areas  for  nonlinear  evaluation  of  the  work  described  herein. 
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Figure  2.6  Simulink  Block  Diagram  Representation  of  the  Engine 
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2.4  Linear  Model  Generation 

In  many  engineering  disciplines  a  physical  process  can  be  modeled  by  a  set  of 
nonlinear  differential  equations  of  the  form: 


x  =  f(x,  u,  t) 
y  =  g(x,  u,  t) 


[2.13] 


where  x  is  the  vector  of  state  variables,  y  is  a  vector  of  outputs,  u  is  a  vector  of  external 
inputs,  t  is  the  scalar  time  variable  and  f  and  g  are  functionals.  These  equations  are 
derived  using  physical  laws  and  are  then  nondimensionalized  so  as  to  minimize  the 
number  of  free  parameters. 

Since  many  processes  are  “weakly”  nonlinear,  they  may  be  approximated  by  a  set 
of  linear,  time  invariant  (LTI)  differential  equations,  called  state  space  models,  given  by: 


x  =  A(a)x  +  B(a)u 
y  =  C(a)x+D(a)u 


where  a  is  a  vector  of  parameters  corresponding  to  the  point  of  linearization  of  [2.13], 
and  A,  B,  C,  D  are  the  system  matrices.  In  the  aircraft  engine  business  it  is  common 
practice  to  use  state  space  models  (Geyser,  1979)  and  (Sellers  and  Daniele,  1975).  These 
models  are  designed  for  the  purpose  of  designing  control  laws  capable  of  safe  operation 
of  the  machinery  while  achieving  high  performance  levels  and  near  optimal  efficiency. 
When  an  engineering  system  can  be  represented  as  in  equation  [2.13]  a  Taylor  series 
expansion  is  formed  at  the  parameter  point  a  as  follows: 
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f(x,  u,  t)  =  f(x(a),  u(a),  t0)  +  |^  (x-x(ct))+|^  (u-u(cx))+...  [2.15] 

OX  x=x(a)  x  OU 

U=lf(a)  X=?I{a) 

ff=TI(a) 

Since  the  system  is  assumed  to  be  in  steady  state  with  the  state  x  =  x(a)  and  control  input 
u  =  u(a),  the  time  derivative  of  the  state  vector  dx/dt  is  zero  and  thus  so  is  the  first  term 
of  the  expansion.  The  series  is  then  truncated  after  the  linear  terms,  resulting  in  the  linear 
system  of  equation  [2. 14].  Here,  the  change  of  variable  x  =  x  -  x(a)  and  u  =  u  -  u(a) 
is  used  for  ease  of  notation. 

When  an  engineering  system  is  not  represented  by  a  set  of  nonlinear  differential 
equations,  but  instead  is  modeled  using  detailed  computer  routines  then  a  perturbation 
technique  must  be  used  in  order  to  extract  linear  models.  A  widely  accepted  technique 
was  developed  (Weinberg,  1975)  and  involves  the  computation  of  the  system  matrices  by 
successively  perturbing  the  model  states  and  measuring  the  responses  within  the 
computer  code.  In  fact  such  a  routine  is  built  into  the  Matlab/Simulink  environment  to 
generate  linear  models.  This  functionality  of  Matlab  was  used  to  create  the  linear  models 
used  in  this  work.  The  nonlinear  model  is  trimmed  at  some  operating  point  {Xq,  U0), 
then  the  linear  model  will  be  a  perturbation  model  of  the  following  form: 

AX  =  A  AX+ B  AU  p16-j 

AY  =  C  AX+ D  AU 
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where  [A,  B,  C,  D]  are  the  partial  derivatives  of  the  state-derivatives  with  respect  to  the 
states  and  the  inputs.  It  is  traditional  to  drop  the  A  from  AX,  AY  and  AU  (note  that 
AU  =  U  -  U0 ,  etc.)  to  obtain  the  linear  state-space  equations  for  the  system. 

In  order,  to  ensure  sound  numerical  properties  of  the  linearized  models,  the  state 
space  system  is  normalized  using  a  scaling  matrix,  R^,  and  a  new  vector  of  state  variables 
z  as  follows: 

x  =  Rxz  [2.17] 

where: 

"rxi  0  O' 

Rx=  0  rx2  0  ,  [2.18] 

.0  0  rx3_ 

and  rxl,  rx2  and  rx3  are  the  A’s  (perturbation  amounts)  used  to  create  the  linear  models.  In 
addition,  the  inputs  and  outputs  are  scaled  by  their  perturbation  values,  defining  a 
normalized  input  vector  u  and  normalized  output  vector  y  defined  by: 

u  =  Ruu,  y  =  Ryy  [2.19] 

where  the  diagonal  matrices  R,,  and  Ry  contain  the  perturbation  amounts  (A’s)  for  the 
inputs  and  outputs,  respectively.  The  state  space  model  then  becomes: 


z  =  R;1  A(a)Rxz  +  Rx'B(a)Ruu 
y  =  R~1C(a)Rxz  +  Ry,D(a)Ruu 


[2.20] 


or  simply: 


z  =  A(a)z+B(a)u 
y  =  C(a)z  +  D(a)u 


[2.21] 


2.5  Plant  Set  Selection 

Plant  sets  considered  for  robust  controller  design  must  satisfy  two  conditions. 
First,  they  must  satisfy  the  path  connectedness  condition.  Secondly,  the  plant  outputs 
must  exhibit  a  reasonable  level  of  response  from  their  associated  inputs  over  the  entire 
region  being  considered.  Evaluation  of  the  sea  level  linear  models  from  idle  to  100% 
corrected  speed  the  thrust  response  to  fuel  flow  input  was  reasonable.  However,  the  stall 
margin  output  to  nozzle  area  input  had  two  distinct  response  characteristics,  this  was 
attributed  to  the  choking  and  unchoking  of  the  nozzle.  This  means  that  the  region  for 
control  design  had  to  be  further  reduced  to  include  only  the  choked  nozzle  region  of 
operation,  which  is  approximately  50%  power  to  100%  power  level. 
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CHAPTER  3 


TRACKING  FILTER  DESIGN 

Model-based  control  will  not  work  well  without  first  having  a  good  model  of  the 
plant  (Adibhatla  and  Gastineau,  1994).  Any  model,  from  linear  models  to  piecewise 
linear  models  to  nonlinear  models,  can  be  used.  However,  to  obtain  the  accuracy  and 
flexibility  required  for  direct  thrust  control,  it  is  useful  to  model  the  physics  of  the  engine 
accurately.  Also,  any  model  that  is  created  will  represent  some  nominal  plant.  Therefore, 
a  means  of  changing  the  model  so  that  it  will  account  for  plant  deterioration  with  age, 
manufacturing  differences,  and  modeling  errors  is  required.  This  is  the  function  of  the 
tracking  filter.  Simply  put,  a  tracking  filter  adjusts  the  model's  inputs  and  parameters  to 
make  it  match  the  actual  plant  outputs.  The  tracking  filter  works  by  forcing  a  selected  set 
of  model  outputs  to  match  the  corresponding  plant  sensors  by  adding  biases  to  selected 
model  inputs  or  parameters.  The  design  is  further  complicated  when  the  plant  can  be 
operated  or  used  in  more  than  one  way.  The  desire  is  that  one  tracking  filter  can  be  found 
that  allows  the  plant  to  operate  in  these  different  modes.  Unfortunately  to  date  such  a 
filter  has  not  yet  been  found. 

3.1  Tracking  Filter  Selection 

The  tracking  filter  is  the  heart  of  any  model-based  control  system.  How  well  it  is 
designed  will  have  a  direct  impact  on  how  well  the  entire  control  system  will  perform.  At 
this  point  it  is  important  for  the  engineer  to  decide  how  the  model-based  control  features 
will  be  used.  Since  every  variable  is  now  available  because  of  the  model  it  appears  that 
the  design  is  simplified.  This  is  not  the  case,  due  to  the  fact  that  the  number  of 
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parameters  that  may  be  updated  or  adjusted  is  much  larger  than  the  number  of  values  that 
can  be  sensed  from  the  actual  plant.  The  fact  is  that  the  difficulty  of  designing  a  control 
has  been  shifted  from  the  closed  loop  controller  design  to  the  design  of  the  tracking  filter. 

The  structure  of  the  tracking  filter  and  model  loop  is  nonlinear.  This  presents 
some  special  design  problems  due  to  the  fact  that  there  is  a  lack  of  design  techniques.  In 
order  to  solve  this  problem  as  suggested  by  Misawa  and  Hedrick  (1989),  the  nonlinear 
model  will  be  linearized  then  the  design  will  proceed  to  find  an  acceptable  design  for 
every  plant  in  the  plant  set.  This  suggests  that  the  use  of  loop  shaping  techniques  would 
be  appropriate  for  the  tracking  filter  design. 

If  the  designer  wants  to  use  a  state-variable  method  and  if  modeling  errors  exist, 
or  the  plant  is  subjected  to  unknown  disturbances,  the  plant  and  model  states  will  be 
different  in  the  steady  state.  As  a  result,  the  corresponding  outputs  will  deviate  from  each 
other.  In  this  circumstance  the  Proportional-plus-integral  (PI)  Observer  is  considered  to 
be  appropriate.  Before  studying  PI  observers,  first  consider  the  development  of  a 
nonlinear  observer  in  an  attempt  to  identify  the  key  issues  that  must  be  addressed  in  any 
design.  Essentially,  it  is  postulated  that  a  dynamic  system  exists  that  has  the  property  that 
the  state  of  the  observer  converges  to  the  state  of  the  observed  process.  The  desire  is  to 
find  this  dynamical  system. 

3.2  Nonlinear  Observers 

An  observer  for  a  plant,  consisting  of  a  dynamic  system 

x  =  f(x,u)  [3.1] 

with  observations  given  by 
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y  =  g(x,u) 


[3.2] 


is  another  dynamic  system,  the  state  of  which  is  denoted  by  x ,  excited  by  the  output  y  of 
the  plant,  having  the  property  that  the  error 

e  =  x-x  [3.3] 


converges  to  zero  in  the  steady  state. 

One  way  of  obtaining  an  observer  is  to  imitate  the  procedure  used  in  a  linear 
system,  namely  to  construct  a  model  of  the  original  system  given  in  equation  [3.1]  and 
force  it  with  the  residual: 


r  =  y-y  =  y-g(x,u)  [3.4] 

The  equation  for  a  proportional  observer  becomes 

k  =  f(x,u)  +  k(y-g(x,u))  [3.5] 

where  k( )  is  a  suitably  chosen  nonlinear  function.  The  differential  equation  for  the  error 
e  can  be  used  to  study  its  behavior.  This  equation  is  given  by 

A 

e  =  x-x 

=  f  (x,  u)  -  f  (x,  u)  -  k(g(x,  u)  -  g(x,  u))  [3.6] 

=  f  (x,  u)  -  f  (x  -  e,  u)  +  k(g(x  -  e,  u)  -  g(x,  u)) 
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Suppose  that  by  the  proper  choice  of  k(  )  the  error  equation  [3.6]  can  be  made 
asymptotically  stable,  so  that  an  equilibrium  state  is  reached  for 

e  =  0 

Then  in  equilibrium,  equation  [3.6]  becomes 

0=  f(x,u)-f(x-e,u)  +  k(g(x-e,u)-g(x,u))  [3.7] 

Since  the  right-hand  side  of  equation  [3.7]  becomes  zero  when  e  =  0,  independent 
of  x  and  u,  it  is  apparent  that  e  =  0  is  an  equilibrium  state  of  equation  [3.6].  This  implies 
that  if  the  k(  )  can  be  chosen  to  achieve  asymptotic  stability,  the  estimation  error  e 
converges  to  zero. 

It  is  very  important  to  notice  that  the  right-hand  side  of  equation  [3.7]  becomes 
zero  independent  of  x  and  u  only  when  the  nonlinear  functions  f(y)  and  g(y)  used  in  the 
observer  are  exactly  the  same  as  in  the  equations  [3.1]  and  [3.2]  that  define  the  plant 
dynamics  and  observations,  respectively.  Any  discrepancy  between  the  corresponding 
functions  will  generally  prevent  the  right-hand  side  of  equation  [3.7]  from  vanishing  and 
hence  will  lead  to  a  steady-state  estimation  error.  Since  the  mathematical  model  of  the 
physical  process  is  always  an  approximation,  in  practice  the  steady-state  estimation  error 
will  generally  not  go  to  zero.  But,  by  careful  modeling,  the  discrepancies  can  usually  be 
minimized  between  the  f  and  g  functions  of  the  true  plant  and  the  model  used  in  the 
observer.  This  can  keep  the  steady-state  estimation  error  acceptably  small. 
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For  the  same  reason,  that  the  model  of  the  plant  and  the  observation  that  is  used  in 
the  observer  must  be  accurate,  it  is  important  that  the  control  signal  u  that  goes  into  the 
plant  is  the  very  same  control  used  in  the  observer.  If  the  control  to  the  plant  is  subject  to 
saturation,  for  example,  then  the  nonlinear  function  that  models  the  saturation  must  be 
included  in  the  observer.  Failure  to  observe  this  precaution  can  cause  difficulties.  This 
problem  is  overcome  by  using  actuator  feedback  values  as  the  inputs  into  the  model, 
thereby,  accounting  for  any  nonlinear  behavior  out  of  the  actuators. 

3.3  Proportional  Plus  Integral  Observer 

Here  we  study  the  proportional-plus-integral  (PI)  observer  that  was  first 
introduced  by  Wojciechowski  (1978)  for  single-input,  single-output  systems  and  further 
developed  for  multivariable  systems  by  Kaczorek  (1979)  and  Shafai  (1985).  Finally, 
Beale  and  Shafai  (1989)  present  an  investigation  of  the  robustness  properties  of  a  PI 
observer  used  in  a  feedback  compensation  configuration.  The  PI  observer  offers  the 
designer  additional  degrees  of  freedom  to  actually  get  the  desired  response. 

Consider  a  linear  time-invariant  system  described  by 

x  =  Ax(t)  +  Bu(t) 
y  =  Cx(t) 

where  the  state  x,  the  input  u  and  the  output  y  are  n-,  m-  and  p-dimensional  vectors, 
respectively.  Without  loss  of  generality  assume  that 
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and  B  has  no  special  form.  Note  that  the  non-zero  entries  are  denoted  by  x  and  A;i,  Q, 
are  V;  x  v;,  v;  x  Vj  ,px  Vj  submatrices  with  V;  being  the  observability  indices  of  the  system. 

It  is  known  that  any  observable  linear  time-invariant  system  of  the  form  given  in 
equation  [3.8]  that  is  not  in  the  observable  canonical  form  can  be  transformed  by  an 
equivalence  relation  to  the  structure  given  by  equation  [3.9],  (Kailath,  1980,  pp.  424- 
437).  The  methodology  for  finding  this  equivalence  relation  is  given  in  Appendix  A. 
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Now,  consider  a  linear  time-invariant  system 


i(t)  =  (A  -  GC)x(t)  +  Gy  (t)  +  Bu(t)  +  Kw(t) 
w(t)  =  y(t)-Cx(t) 


where  x(t)  and  w(t)  are  n-  and  p-  dimensional  vectors,  respectively.  G  and  K  are 
matrices  of  the  proportional-plus-integral  observer  that  must  be  determined. 

Definition  3.1  The  system  described  by  equation  [3.11]  is  said  to  be  a  full-order  PI 
observer  for  the  system  given  in  equation  [3.8]  if  and  only  if: 


lime(t)  =  0 

t-»cO 

limw(t)  =  0 


[3.12] 


where  e(t)  =  x(t)  -  x(t)  represents  the  observer  error. 

Theorem  3.1  (Beale  and  Shafai,  1989)  The  system  given  in  equation  [3.11]  is  a  full- 
order  PI  observer  for  the  system  given  in  equation  [3.8]  if  and  only  if  all 
the  eigenvalues  of  the  matrix: 


R  = 


A-GC  K 
-C  0 


[3.13] 


have  negative  real  parts. 

Proof.  From  equations  [3.8]  and  [3. 1 1]  the  dynamics  of  the  observer  error  becomes: 

e(t)  =  (A  -  GC)e(t)  +  Kw(t)  [3.14] 


38 


which  yields 


'e(t)' 

'A-GC 

K' 

e(t) " 

w(t). 

-C 

0_ 

w(t)_ 

[3.15] 


and  it  follows  immediately  that  x(t)  is  an  estimate  of  x(t)  if  and  only  if  Re  Xq[R]  <  0 
for  q  =  l,...,  (n  +  p). 

Thus  the  problem  reduces  to  one  of  choosing  the  matrices  G  and  K  such  that  the 
system  given  in  equation  [3.11]  will  be  a  PI  observer  for  the  system  given  in  equation 
[3.8].  A  methodology  for  selecting  the  matrices  G  and  K  for  both  single-output  and 
multiple-output  systems  is  known.  These  methodologies  were  presented  by  Beale  and 
Shafai  (1989). 

3.4  Proportional-Plus-Integral  Observer  Design  for  Single-Output  Systems 

In  this  case  p  =  1  and  the  system  matrices  have  the  form 
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Let  the  unknown  PI  observer  matrices  G  and  K  be  written  as 
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[3.17] 


Then  a  simple  algebraic  manipulation  shows  that 


M„-(A-GC) 

c 


-K 


g0-l) 


[3.18] 


and  by  comparison  of  its  coefficients  with  the  ones  specified  by  the  desired  characteristic 
equation,  i.e. 


A„M=nVx>)  p-19] 

i=l 

this  results  in  a  system  of  (n  +  1)  equations,  which  have  to  be  satisfied  by  2n  elements  of 
G  and  K.  Therefore,  (n  -  1)  elements  of  G  and  K  can  be  chosen  arbitrarily  (except  g(n.,) 
and  lq,)  and  the  remaining  elements  can  be  evaluated  from  the  system  of  (n  +  1) 
equations. 
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3.5  Proportional-Plus-Integral  Observer  Design  for  Multiple-Output  Systems 

In  the  general  case  of  p  outputs,  we  decompose  the  system  given  in  equation  [3.8] 
into  p  subsystems  each  having  dimension  v,  as  specified  by  equations  [3.9]  and  [3.10]. 
Let  the  state  vector  of  each  subsystem  be  represented  by 


X;  = 


«<M)  +  1 


C0-2) 


X„ 


+  1 


,  i  —  1,  p 


[3.20] 


where 


[3.21] 


Then  the  ith  subsystem  can  be  written  as 


F 

Xj  =Aiixi+^Aijxj+Biu 


j=i 


[3.22] 


or  equivalently  as 


F 

xi=Ajixi+Xaioxa.+Biu 


j=i 


[3.23] 
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where  An,  Ay,  and  Bj  are  the  submatrices  defined  by  equations  [3.9]  and  [3.10],  and  a,aj  is 
the  CTjth  column  of  the  matrix  Ay.  Also  let  Cc  be  the  pxp  lower  triangular  matrix  formed 
from  C  by  eliminating  all  by  the  p  non-trivial  columns  of  C  so  that  we  have 


1  0  -  0 

x  1  •••  0 


x  x  •••  1 


[3.24] 


Using  [3.8]  we  generate  a  new  output  equation  as 


y  =  c;‘y = 


x 


<*2 


X„ 


[3.25] 


so  that  the  XCTi,  i  =  1,  ...,  p,  are  available.  Thus,  it  is  apparent  that  the  ith  subsystem 
described  by  equation  [3.23]  is  a  vrdimensional  single-output  system  in  observable 
companion  form  driven  by  the  directly  measurable  system  signals  u(t)  and  xc.  (j  ^  i). 
Therefore,  the  procedure,  which  was  given  for  single-output  systems,  can  now  be 
employed  for  each  of  the  p  subsystems. 


3.6  Linear  Versus  Nonlinear  Tracking  Filter 

The  design  of  a  nonlinear  tracking  filter  is  possible.  It  would  most  likely  be  based 
upon  neural  networks  or  fuzzy  logic  techniques.  These  techniques  have  recently  gained 
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favor  within  the  control  community.  But,  the  amount  of  data  necessary  to  design  a  neural 
network  was  considered  prohibitively  large  for  the  task  that  needed  to  be  accomplished. 
For  this  work  a  linear  tracking  filter  was  deemed  sufficient  for  the  following  reasons. 
First,  the  dynamics  of  the  engine  change  slowly  as  it  is  flown  throughout  the  flight 
envelope.  Secondly,  the  error  between  the  actual  plant  and  the  model  will  in  general  be 
small.  This  is  due  to  the  amount  of  time  developing  the  model  during  the  engine 
development  process.  This  small  error  means  that  any  correction  to  the  model  will  adjust 
the  model  outputs  in  a  linear  fashion.  Finally,  linear  observer  techniques  are  better 
understood  and  offer  a  satisfactory  solution  for  this  problem. 

If  the  plant  had  been  highly  nonlinear  or  the  dynamic  behavior  had  changed 
rapidly  then  the  only  choice  may  have  been  to  design  a  nonlinear  tracking  filter.  Also,  if 
the  model  contains  first  principle  physics  of  the  plant  and  is  not  necessarily  corrected  or 
adjusted  based  upon  development  testing  results  then  nonlinear  techniques  may  eliminate 
or  reduce  the  development  time  required  for  the  model. 
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CHAPTER  4 


QUANTITATIVE  NYQUIST  ARRAY  USING 
FORWARD  PATH  DECOUPLING 

The  idea  of  interaction  between  control  loops  is  of  paramount  concern  when 
developing  control  systems  for  multiple-input,  multiple-output  systems.  For  any  practical 
design  to  be  decoupled  or  to  have  very  little  interaction  between  loops  is  highly  desirable 
(Suh,  1990).  In  this  chapter,  a  methodology  for  reducing  loop  interaction  is  presented 
and  then  tied  to  a  quantitative  design  for  a  diagonal  feedback  controller. 

After  the  tracking  filter  has  been  successfully  designed  and  tested  to  meet 
specifications  one  can  now  consider  the  n-input,  n-output  feedback  control  system  shown 
in  Figure  4.1.  In  a  model-based  control  structure,  the  dynamics  of  the  tracking  filter  and 
model  are  now  included  in  the  definition  of  the  plant  P.  The  feedback  values  are  now 
model  computed  values. 


Figure  4.1  Feedback  Structure 
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The  transfer  function  for  the  system  shown  above  is  given  by: 


ty/r  =  (i+pkg)'‘pkgf 
=  (p+kg)"’kgf 

1 

where  P  =  P  1  and  has  the  elements  p;j  =  — . 

*lij 


[4.1] 


The  problem  that  is  being  considered  here  is  the  quantitative  feedback  design  problem 
specified  by  Horowitz  (1979)  of  designing  a  decoupling,  precompensating  matrix,  K, 
with  diagonal  elements  equal  to  unity  (kH  =  1),  a  diagonal  feedback  controller,  G,  and  a 
prefilter,  F,  to  achieve  certain  tracking  specifications  given  by: 


Aij(£o)<TY/R(jco)<Bij(»)  i,j  =  l,  2,  ...n  [4.2] 


Ajj(co)  and  B^co)  are  the  client’s  desired  response  characteristics,  for  all  plants,  PeP, 
where  due  to  uncertainty,  P  is  the  set  of  all  possible  plants.  All  matrices  are  nxn. 

The  structure  of  Figure  4.1  was  selected  over  the  structure  used  by  Nordgren, 
etal.  (1994)  because  of  the  apparent  ease  of  implementation.  The  approach  used  by 
Nordgren  emphasized  the  ease  of  design  over  implementation.  However,  it  is  important 
to  note  that  the  level  of  performance  achievable  by  the  two  configurations  is  identical. 
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4.1  Preliminary  Mathematics 

The  following  is  a  collection  of  definitions  and  elementary  results  that  will  be 
used  for  the  work  contained  herein.  It  is  based  on  the  theory  of  non-negative  matrices 
(Berman  and  Plemmons,  1979). 

Definition  4.1:  A  matrix,  M  e  Rnxn  is  an  M-matrix  if  the  diagonal  elements  of  M  are 
positive,  the  off  diagonal  elements  of  M  are  non-positive,  and  the 
principal  minors  of  M  are  non-negative. 

Definition  4.2:  Given  a  matrix  Z  e  Cnxn  the  comparison  matrix  of  Z,  M(Z),  has  elements, 


Definition  4.3: 


(diagonal  elements) 

i 56  j 

A  matrix,  Z  eCnxn  is  irreducible  if  there  does  not  exist  a  permutation, 
P  eCnxn  such  that 


M(Z)jj  =  • 


I 

1- 

zn 

Zj; 

L 

u 

Definition  4.4: 

Definition  4.5: 
Definition  4.6: 


PZP-1  = 


7  1  7 

_  n_}_  1 

0  1  7 

V  |  ^22 , 


,  with  square  submatrices,  Zn  and  Z 


22* 


An  irreducible  matrix,  Z  e  Cnxn  is  an  H-matrix  (Hadamard  matrix)  if 
M(Z)  is  an  M-matrix. 

An  H-matrix,  Z  e  Cnxn  is  a  BH-matrix  if  Z'1  is  also  an  H-matrix. 

A  matrix,  Z  admits  a  diagonal  regular  splitting  if  it  can  be  written  as 
Z  =  D  +  C ,  with  D  =  diag{Z}  non-singular. 
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Definition  4. 7:  The  Perron  root  of  a  matrix  M,  is  the  largest  eigenvalue  of  the  matrix 
formed  from  the  absolute  values  of  the  elements  of  M. 

?ip(M)  =  max{^(M+)}  =  p(M+) 

p(M+)  is  the  spectral  radius  of  M+.  The  absolute  values  are  taken 
element-wise,  (M+)..  =|My|. 

The  Perron  root  is  analytic  and  monotonic  with  respect  to  the  elements 
of  M+.  The  Perron  root  is  an  upper  bound  on  the  eigenvalues,  A,(M),  and 
therefore  the  spectral  radius,  p(M),  of  a  matrix: 

|X.(M)|  <  p(M)  <  A,p(M) 

The  spectral  radius  is  a  lower  bound  on  all  induced  norms. 

Definition  4.8 :  Given  that  Z  eCnxn  admits  the  diagonal  splitting, 

Z  =  D  +  C  =  (l  +  CD'1  )D  =  (I  +  M)D 

M  is  the  interaction  matrix  and  the  interaction  index  of  Z  is  the  Perron 
root  of  the  interaction  matrix, 

y(z)-Mm) 
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The  utility  of  the  interaction  index  for  feedback  design  is  seen  by  considering  the  inverse 


of  Z,  (Z  1  =  Z).  If  Z  is  an  H-matrix,  M  has  spectral  radius  less  than  unity 
(p{M)<Xp(M)  <  1,  with  y(Z)  =  Zp(M))  and  (l  +  M)_1  =j^(-M)k  converges. 

'  ‘  k=0 

Also, 


K(z ) 


Xp(Z)  =  »,(D1UM)> [43] 


Definition  4.9 :  An  irreducible  matrix  Z  <=  Cnxn  is  almost  decoupled  if  the 

max|A,p  (Z),  Xp  (z)|  <  0.5.  (Z  almost  decoupled  =>  Z  is  a  BH  matrix  - 
see  below.) 

Theorem  4. 1 :  An  irreducible  matrix,  Z  €  Cnxn  is  a  BH-matrix  if 

min|xp(Z),  Ap  (zj|  <  0.5.  Notice  that  the  concept  of  almost  decoupling 

is  stronger  than  the  requirements  for  Z  to  be  a  BH-matrix. 

Proof.  See  (Nwokah,  et  al,  1995)  or  appendix  B. 


Mz) 


Theorem  4.2 :  An  H-matrix  Z  eCnxn  is  almost  decoupled,  if  X.p(z)  <  - — v  <  0.5; 


requiring  ?ip(Z)<l/3. 

Proof.  This  proof  follows  from  equation  [4.3]  and  definition  4.9. 
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Theorem  4.3:  Given  P  eP,  square  and  invertible  with  no  hidden  unstable  modes,  the 
system, 

Ty/r  =  (i+pkg)~'pkgf 

=  (P0  +  G)'1  (i  +  (P0  +  (K  -  I)G)(PD  +  G)"  )"  KGF 
=  (pd+G)_,(I  +  M)''KGF  [4.4] 

is  internally  stable  if: 

i)  G,  F,  K  are  designed  to  be  stable. 

ii)  Xp(M)<l. 

iii)  Each  (1  +  giqH )  is  designed  to  have  no  zeros  in  the  right  hand  plane 
For  proof  see  appendix  B. 

4.2  Insights  of  Interaction  in  Multiple-Input,  Multiple-Output  Systems 

The  idea  of  being  decoupled  or  having  very  little  interaction  between  loops  is 
highly  desirable  in  any  practical  design  (Suh,  1990).  The  goals  of  decoupling  can  be 
illustrated  as  follows.  Consider  the  return  difference  matrix  and  its  regular  splitting  into  a 
diagonal  matrix  D  and  an  off-diagonal  matrix  C: 

(I  +  PKG)  =  D  +  C  =  (i  +  CD_I)D  [4.5] 
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The  goal  of  decoupling  is  to  force  the  spectral  radius  (p)  of  the  modulus  of  the  product 
CD'1  to  be  less  than  one.  In  this  way,  the  diagonal  elements  of  the  return  difference 
matrix  solely  determine  the  closed  loop  stability  of  the  system. 

Methods  such  as  the  Relative  Gain  Array  (RGA)  are  used  to  measure  the 
interaction  between  inputs  and  outputs.  It  is  shown  in  Section  4.1  that  the  Perron  root  can 
be  a  useful  measure  of  the  degree  of  interaction  in  a  feedback  design.  It  is  important  to 
point  out  that  the  interaction  index  is  a  measure  of  generalized  diagonal  dominance  rather 
than  strict  diagonal  dominance  of  a  matrix. 

4.3  Precompensation  Matrix  Design 

Consider  equation  [4.1]  and  separate  it  so  that  the  diagonal  and  off-diagonal 
components  are  identified. 

ty/r  =(p+kg)_1kgf 

=  (PD  +  P0  +  (K  -  I)G  +  G)“‘  KGF  , 

=  (pd  +  g+p0+(k-i)g)_,kgf 

=  (PD  +  G)'1  (i  +  (P0  +  (K  -  I)G)(Pd  +  G)"’ )  ’  KGF  [4.6] 

Notice  that  equation  [4.6]  is  in  the  form  of  equation  [4.4]  where: 

C  =  P0  +(K-I)G 
D  =  (Pd+G) 
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therefore; 


M  =  CD'1  =  (P0  +  (K-I)G)(PD  +  a)-1' 


[4.7] 


Now  for  a  given  plant  set  it  is  necessary  to  determine  how  to  design  K  so  that  the 
interaction  index,  y  =  >.p(M) ,  is  less  than  unity.  First,  define  the  diagonal  sensitivity  (S) 
and  diagonal  complementary  sensitivity  (T)  as: 


T  =  G(pD+G)-',tB=T|^-  =  Tfi- 


[4.8] 

[4.9] 


The  perron  root  is  found  by  solving  the  maximum  eigenvalue-eigenvector  problem  of 
equation  [4.7]  which  is  as  follows. 


(p„+(k-i)g)(pd+g)' 
|p0(pd+g)"‘+(k-i)g(pd+g)' 
p„((i+gp;%)%(k-i)g(pd+g) 
IPoP^I  +  GPJ1)"'  +(K-I)G(Pd  +g) 


x  =  A.px 
x  =  Xpx 
x  =  Xpx 
x  =  Xpx 


[4.10] 


Substituting  equations  [4.8]  and  [4.9]  into  equation  [4.10]  the  following  result  is 
obtained. 
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|p0pd-'s+(k-i)t|x=|n|x=».px 


[4.11] 


with  the  elements  of  N  given  by: 


nu  = 


<lji  1 


qsl  +  <j+k*l  +  <J 


A 


_kU 
V'iij  J 


1 


+  kjj 

\  +  t,  ,J 


[4.12] 


If  only  the  specification  for  the  magnitude  of  the  sensitivity  of  individual  loops  is 


known,  i.e. 


following. 


1  +  t  ! 


<Sj(ni),  we  may  apply  the  Schwartz  inequality  to  obtain  the 


ny  * 


TJJ 


qy 


-kij 


Isj+ky 


qy 


-k: 


+  k, 


[4.13] 


To  guarantee  reduction  of  the  interaction  index  using  the  precompensating  matrix, 
K,  which  has  only  off-diagonal  elements  for  design,  may  be  reduced  by  designing 


^jj 


-ky 


^jj 


qy 


or 


max 
P  eP 


1-k  — I 

ijqJ 


<1. 


[4.14] 
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This  can  be  achieved  over  frequency  using  a  QFT  type  design.  However,  from  equation 
[4.12]  the  following  ideas  need  to  be  considered  when  designing  the  matrix  K.  First, 

when  |*j| « 1  we  find  that  jn^  « jq^/q^j  so  that  the  interaction  index  is  the  same  as  that 

of  the  open  loop  plant.  Secondly,  when  |^|»1  then  we  find  that  |nj*|k8|.  This 
motivates  choosing  the  off-diagonal  elements  of  K  as  small  as  possible.  Finally,  the  third 
consideration  is  that  when  |^| » 1 ,  especially  around  the  gain  and  phase  cross  over 


frequencies  of  where  typically  a  margin  such  as 


1 

1  +  t} 


*P. 


P  >  0  dB  is  specified, 


then  selecting  k-  »  —  will  keep  individual  elements  of  |N|  small. 

% 

4.4  Diagonal  Loop  Controller  Design 

There  exist  several  methodologies  for  a  multivariable  QFT  design.  The  technique 
used  here  makes  use  of  the  plant  inverse.  Now  consider  equation  [4.1],  after  the 
precompensating  matrix,  K,  has  been  designed  it  is  easier  to  rework  the  equation  into  the 
following  form. 


ty/r  =(p+kg)“'kgf 
=  (kp+g)"’gf 


[4.15] 


Then  defining 


one  has  the  standard  QFT  design  problem  given  by: 
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[4.16] 


Tv,e=(p'+g)"’gf 

Using  simple  algebraic  manipulations  equation  [4.16]  is  rearranged  into  the  following 
form,  then  specific  equations  are  derived  for  each  ty  in  the  matrix  of  transfer  functions. 

(P'+G)Ty,r=GF  [4.17] 

In  order  to  illustrate,  consider  the  2x2  system  with  G  diagonal  and  F  fully 
populated.  Expanding  equation  [4. 1 7]  in  detail  one  gets: 

IT—  — 

qli  q« 

V^_q.21  ^22 

Equation  [4.18]  contains  four  multiple-input,  single-output  equations  given  by: 


All  four  equations  can  be  reduced  to  the  following 
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[4.20] 


l  Aju _ L_Snt 

‘“"l+O+Aq*  21 

_  4^21 _ 1  4  22  j. 

71  1  +  4  1+^2  q2i 


_  4jl2 _ 1  Qll 

1  +  4  1  +  4  q12 

_  44 1 _ 1  0.22 

1  +  4  1  +  4  q21 


where  4  =  g^ . 

Or  generically  this  can  be  written  as: 


f-l 

1  ( 

1  +  £: 

1  +  4 

q« 


forj  =  l...m. 


[4.21] 


Now  specifications  on  the  individual  loop  sensitivities  can  be  established  to  reduce  the 
interaction  index  to  required  levels.  These  specifications  would  be  in  addition  to  the 
uncertainty  reduction  and  disturbance  rejection  specifications  that  are  used  in  standard 
QFT  designs. 

4.5  Prefilter  Design 

In  order  to  achieve  enhanced  or  faster  response,  a  prefilter  is  placed  at  the  input  to 
each  channel.  The  specifications  for  most  QFT  designs  are  generally  non-interacting. 
The  feedback  controller  is  designed  to  achieve  specified  levels  of  uncertainty  reduction 
for  the  diagonal  closed  loop  elements.  Then  the  prefilter  would  only  have  diagonal 
elements  designed  such  that  the  closed  loop  response  would  lie  within  the  client  s 
specifications.  If  the  closed  loop  system  is  lower  triangular  or  almost  lower  triangular 
then  we  may  be  able  to  further  reduce  the  open  loop  interaction  index  by  designing  off 
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diagonal  elements  in  the  prefilter.  Consider  the  following  example  to  illustrate  the 
methodology  for  designing  the  prefilter. 

Suppose  T  =  (I  +  L)'1  L  =  8,2  ,  with  |sijj«l.  Choose  F  =  11  and  design 

l_f  21  f  22  J  L^21  M2  _ 

Ty/R  ~  TF  in  the  following  order. 


1)  First  design  the  diagonal  elements  fH  to  achieve  the  client’s  tracking  specifications  on 
(Tm)^  i  =  1, 2,  3. 


2)  For  a  reduction  in  interaction, 


j(TY/R  )2]  |  —  1*21^11  +  ^22^21 1  —  |*2l|’  0r  ^+.  f  ^21  -  r- 

MIMl  Ml 


[4.22] 


Note,  in  the  worst  case,  the  same  interaction  index  can  be  obtained  by  choosing  fi}  =  0. 
This  design  can  be  accomplished  using  techniques  from  Quantitative  Feedback  Theory. 
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CHAPTER  5 


DESIGN  EXAMPLE 


5.1  Engine  Model 

The  equations  used  in  Chapter  2  were  used  to  develop  a  nonlinear  turbojet  engine 
simulation  model  for  the  purpose  of  designing  a  single  multivariable  controller  for  an 
engine  at  sea  level  static  conditions.  This  model  was  then  linearized,  resulting  in  thirty- 
one  plant  models  or  sets  of  state  space  matrices.  Refer  to  appendix  D  for  this  set  of 
normalized  plant  models.  Each  of  these  linearized  models  was  extracted  from  the 
nonlinear  model  using  a  perturbation  technique.  Each  model  consists  of  two  inputs,  three 
states,  and  either  two  outputs  if  designing  the  controller  and  precompensating  matrix  or 
three  outputs  if  designing  the  tracking  filter.  This  is  summarized  in  Table  5.1.  Also 
included  to  improve  the  realism  of  the  problem  were  models  of  the  transport,  combustion, 
and  sensor  time  delays,  as  well  as  actuator  models  for  nozzle  exit  area  control  and  fuel 
flow  rates.  A  schematic  diagram  of  these  details  is  shown  in  Figure  5.1. 

In  order  to  obtain  a  successful  robust  controller  design,  the  path  connectedness 
condition  must  be  satisfied,  and  all  the  outputs  must  exhibit  a  reasonable  amount  of 
response  due  to  their  associated  inputs  (Nordgren,  1994).  The  linear  models  used  for  this 
work  were  trimmed  about  values  of  percent  corrected  speed  that  ranged  from  70%  (idle 
speed)  to  100%  (maximum  rotation  speed).  A  study  of  path  connectedness  was 
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Table  5.1  Engine  Model  Variables 


Outputs 

oTclI&S 

Tracking  Filter 

Control 

Speed 

Speed 

Thrust 

Combustor 

Pressure 

Compressor  Exit 
Pressure 

Stall  Margin 

Tailpipe 

Pressure 

Turbine  Exit 
Pressure 

Figure  5.1  Feedback  Configuration  Showing  Complete  Modeling  Details  After  the 
Tracking  Filter  Design  is  Complete. 

performed  and  plant  number  one  through  number  twenty-nine  were  found  to  be  path 
connected.  Closer  evaluation  of  plants  thirty  and  thirty-one  revealed  that  the  perturbation 
levels  used  to  generate  the  linear  models  actually  resulted  in  the  model  interpolating  a 
value  for  the  compressor  that  was  off  the  map.  Therefore,  there  is  not  a  necessarily 
smooth  connected  path  set  when  we  move  off  the  maps,  this  resulted  in  the  elimination  of 
these  plants  from  the  acceptable  set.  This  path  connectedness  study  made  use  of  the  ratio 
of  the  determinants  of  the  plant  nominal  model  and  the  specific  plant  of  interest.  If  this 
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plant31 


Figure  5.3  Ratio  of  Determinants  for  a  Nonpath  Connected  Plant 

To  further  improve  the  overall  design  the  acceptable  plant  set  is  reduced  by 

investigating  the  open  loop  plant  responses  when  perturbed  using  a  unit  step.  Figure  5.4 

shows  the  Thrust  and  Stall  margin  responses  to  individual  inputs  of  fuel  flow  and  nozzle 

area.  Note  that  for  the  thrust  response  to  fuel  flow  input  that  characteristics  tend  to  group 

into  two  groups,  one  with  a  fast  time  response  and  the  other  with  a  slower  response.  This 

characteristic  continues  for  the  stall  margin  response  to  a  change  in  nozzle  area. 

However,  it  is  not  as  pronounced.  This  means  that  there  is  more  uncertainty  in  the 

response  to  the  input.  Because  of  this  the  thrust  response  is  used  to  identify  the  subset  of 

acceptable  plants  to  use  and  the  stall  margin  response  will  be  worked  with.  This  choice 
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could  be  different  if  engine  performance  was  not  of  primary  interest.  Figure  5.5  shows 
the  open  loop  responses  for  the  acceptable  subset  of  plants.  These  are  the  plants  that  the 
subsequent  design  was  done  on.  Note  that  the  additional  plants  are  not  forgotten.  After 
the  design  is  complete  these  plants  will  be  evaluated  to  determine  how  well  they  meet  the 
specification.  Then  if  need  be,  the  specification  can  be  changed  or  a  new  design  can  be 
performed  for  this  set  of  plants.  This  is  the  traditional  engineering  problem. 


Delta  Fuel  Flow  Delta  Nozzle  Area 


Delta  Fuel  Flow  Delta  Nozzle  Area 


Figure  5.5  Open  Loop  Responses  for  the  Family  of  14  Plants 


5.2  Control  Objective  and  Specifications 

The  control  objective  is  to  regulate  the  two  engine  outputs  (Thrust  and  Stall 
margin)  that  are  unmeasurable  using  the  two  independent  variables  (Fuel  Flow  and 
Nozzle  area)  to  achieve  some  overall  response  specification.  For  this  problem  there  are 
different  response  characteristics  for  each  loop.  A  time  domain  envelope  described  by 
both  an  upper  and  lower  second  order  response  models  provides  the  regions  of  acceptable 
design.  These  models  are  given  by: 
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[5.1] 


A  nominal  plant  model  was  chosen  arbitrarily  to  be  plant  number  sixteen  out  of 
the  family  of  twenty-nine  plants.  This  corresponds  to  the  engine  being  at  sea  level  static 
and  a  percent  corrected  rotor  speed  of  85%.  The  normalized  nominal  models,  without  the 
computational  delays,  actuator  and  sensor  dynamics  and  transport  delays  is  given  in  state- 
space  form  for  the  control  and  tracking  filter  designs  in  equations  [5.2]  and  [5.3] 
respectively.  These  equations  differ  only  by  the  output  equations. 
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5.3  Tracking  Filter  Design 

The  tracking  filter  was  designed  using  the  methodology  described  in  Chapter  3. 
A  tracking  filter  was  sought  that  converged  the  model  outputs  to  the  plant  outputs  quickly 
while  at  the  same  time  providing  good  thrust  and  stall  margin  computations. 
Unfortunately,  there  exists  no  specification  for  the  tracking  filter  response.  Therefore,  a 
specification  was  established  that  the  eigenvalue  (or  pole  locations)  should  be  located  so 
that  a  second  order  response  achieved  steady-state  approximately  twice  as  fast  as  the 
fastest  responding  plant  within  the  plant  set.  This  implies  that  the  design  should  be  done 
with  the  plant  that  has  the  slowest  response.  While  all  the  plants  had  approximately  the 
same  transient  response,  it  was  noticed  that  the  steady-state  error  was  greater  at  lower 
speeds.  Therefore,  a  lower  speed  condition  was  selected  in  order  to  guarantee  that  the 
proportional  plus  integral  observer  offered  excellent  results  at  all  conditions.  Plots  of 
engine  and  engine  model  outputs  before  tracking,  are  shown  in  Figure  5.6. 

After  the  successful  design  of  the  Proportional  plus  Integral  (PI)  Observer  the  PI 
elements  for  each  loop  were  found  to  be. 


PI1 
PI2 
PI  3 
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31s  +  45 


[5.4] 
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In  Figure  5.7  a  comparison  of  the  outputs  is  given  with  the  tracking  filter  switched 
on.  These  results  indicate  that  the  computation  of  thrust  and  stall  margin  will  be  good  at 
all  conditions. 
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Figure  5.6  Outputs  Given  the  Tracking  Filter  is  Off  (-  engine,  —  model) 


Tracking  Filter  is  On  (-  engine,  --  model) 


An  important  fact  to  note  is  that  all  outputs  were  not  required  to  be  tracked  with 
zero  steady-state  error.  In  fact  some  parameters  will  not  track  in  order  to  guarantee 
tracking  of  the  outputs  of  interest.  This  has  significant  implications  if  the  model  is  used 
for  purposes  other  than  achieving  engine  performance.  Such  an  example  is  engine  health 
monitoring  where  the  tracking  of  component  life  is  important.  In  this  case  accurate 
knowledge  of  the  temperature  levels  that  the  components  have  been  exposed  to  have  a 
dramatic  impact  on  the  life  remaining  in  the  component.  Therefore,  temperatures  would 
not  be  allowed  to  have  as  large  an  error  as  shown  in  Figure  5.7. 

In  Figure  5.8  a  step  in  the  demanded  fuel  flow  was  introduced  to  determine  how 
well  the  model  would  then  track  the  plant  after  the  steady  state  has  been  achieved.  This 
plot  demonstrates  that  once  the  system  is  matched  then  it  is  consistent  given  a  disturbance 
or  movement  to  another  operating  condition. 

5.4  Forward  Loop  Design 

After  the  tracking  filter  has  been  designed,  the  control  design  can  begin.  Control 

system  design  becomes  somewhat  trivial  when  all  parameters  are  available.  However,  in 

this  example  where  thrust  and  stall  margin  are  being  controlled  simultaneously  there  is 

significant  amounts  of  interaction  between  the  inputs  and  the  outputs.  Thereby 

complicating  the  design.  The  control  structure  of  Chapter  4  is  chosen,  which  includes  a 

decoupling  precompensating  system  and  a  diagonal  closed  loop  controller.  Finally,  a 

prefilter  is  necessary  to  shape  the  output  responses  such  that  the  time  domain  response 

specifications  are  satisfied  to  the  degree  feasible.  In  addition,  the  dynamics  or  time  delay 

associated  with  the  updating  of  the  model,  and  the  time  required  for  the  model  to  provide 

an  output  for  feedback  must  be  accounted  for.  Initial  investigations  indicate  that  due  to 
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Tracking  Filter  is  On  with  A  Step  in  One  Input  (-  engine,  --  model) 


time  delays  in  the  sensors  and  actuators  the  maximum  achievable  loop  gain  cross-over 
frequencies,  over  all  the  loops  and  over  the  uncertainty  range  of  the  plant,  and  at  which 
we  expect  the  interaction  will  be  critical,  will  be  in  the  range  of  1 .0  rad/s  to  6.0  rad/s. 

5.4.1  Precompensator  Design.  The  role  of  the  precompensating  system  is  to  reduce  the 
interaction  index  of  the  plant  as  much  as  possible.  The  interaction  index  for  the  plant  is 
given  in  Figure  5.9.  At  low  frequencies  the  interaction  index  is  nearly  less  than  a  third 
indicating  that  the  system  is  almost  decoupled.  However,  for  the  bandwidth  of  concern 
some  of  the  plants  have  interaction  index  greater  than  one,  therefore,  stability  cannot  be 
guaranteed  for  these  plants. 


Figure  5.9  Interaction  Index  for  Open  Loop  Plants 
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The  precompensating  system  has  unity  on  the  diagonal  and  the  elements  to  be 
designed  on  the  off-diagonal.  These  elements  were  originally  designed  to  be  static 
compensators.  But,  due  to  a  push-pop  effect  it  was  found  that  static  compensators 
actually  made  the  interaction  index  worse  at  low  frequencies.  This  was  unacceptable, 
therefore,  a  dynamic  compensator  was  designed.  The  off-diagonal  elements  of  the  pre¬ 
compensator  are  given  as: 
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Figure  5.10  shows  the  resulting  plot  of  interaction  index.  Notice  that  for  all 
frequencies  within  the  bandwidth,  the  interaction  index  is  less  than  one,  thereby  stability 
is  guaranteed.  Also,  through  the  known  system  bandwidth  the  interaction  index  was 
approximately  one-third  or  less,  so  that  the  almost  decoupling  condition  has  been 
achieved.  This  greatly  simplifies  the  controller  design  covered  in  the  next  section. 

5.4.2  Diagonal  Controller  Design.  The  methodology  and  equations  presented  in 
Chapter  4  were  used  as  the  basis  for  designing  the  diagonal  controllers.  This 
methodology  results  in  the  least  conservative  design  in  that  they  utilize  all  available 
parametric  uncertainty  information  available.  For  this  design,  it  did  not  matter  which 
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loop  you  started  with.  However,  loop  one  was  selected  due  to  the  tighter  behavior  of  the 
step  responses  shown  in  Figure  5.5. 


Figure  5.10  Open  Loop  Interaction  Index  for  Plant  and  Precompensator 
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5.4.3  Prefilter  Design.  In  order  to  achieve  a  faster  cleaner  response,  a  prefilter  is  placed 
at  the  input  to  each  channel.  The  prefilter  is  designed  to  force  the  individual  channel 
closed  loop  frequency  responses  to  lie  between  the  upper  and  lower  step  response 
bounds.  Figures  5.11  and  5.12  show  the  frequency  responses  of  the  individual  channels 
plotted  against  their  respective  upper  and  lower  step  response  specifications.  These 
prefilters  are  given  as  follows: 
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5.5  Design  Analysis 

The  resulting  closed  loop  time  domain  step  responses  for  the  feedback  system  are 
shown  in  Figures  5.13  and  5.14.  The  responses  are  presented  in  columns  where  column  1 
shows  the  two  outputs  for  a  step  input  on  channel  one.  Column  two  shows  the 
corresponding  outputs  for  a  step  input  on  channel  two.  Figure  5.13  shows  the  output 
responses  for  the  feedback  system  without  the  prefilter,  and  Figure  5.14  shows  the 
corresponding  outputs  with  the  prefilter.  The  first  channel  shows  a  settling  time  on  the 
order  of  two  seconds,  and  the  second  channel  settling  time  is  on  the  order  of  five  seconds. 
Channel  one  appears  to  have  the  most  uncertainty.  This  is  attributed  to  the  significantly 
larger  time  delay  present  in  this  loop. 
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Figure  5.12  Frequency  Response  for  Stall  Margin  to  Nozzle  Area  Channel 
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CHAPTER  6 


CONTRIBUTIONS,  CONCLUSIONS  AND  RECOMMENDATIONS 
6.1  Summary  of  Contributions 

The  model-based  control  concept  is  a  significant  departure  from  the  typical 
aircraft  engine  control  system  design.  Until  recently,  the  computing  power  available  to 
the  engine  control  designer  has  been  one  of  the  main  limiting  factors.  By  utilizing  the 
improvements  in  today’s  computing  equipment,  the  functionality  of  the  engine  control 
can  be  significantly  increased.  The  model-based  control  concept  capitalizes  on  this 
increase  in  computing  power  so  that  the  designer  can  control  the  engine  in  ways  never 
before  possible.  The  work  presented  herein  contributes  in  several  ways  to  the  overall 
knowledge  of  engine  multivariable  control  system  design.  The  ideas  presented  are  shown 
in  Chapter  5  to  be  very  effective  when  applied  to  a  turbojet  design  problem. 

Model-based  control  (MBC)  can  be  categorized  as  an  adaptive  control  technique. 
However,  model-based  control  differs  from  classical  adaptive  control,  one  of  the  main 
techniques  is  model-reference  adaptive  control  (MRAC),  in  at  least  two  key  areas.  The 
first  difference  is  with  the  controllers.  In  an  MRAC  system  the  parameters  of  the 
controller  are  adjusted  to  obtain  a  desired  response  while  in  a  MBC  system  the  controller 
parameters  are  fixed.  The  second  difference  and  probably  the  most  significant  difference 
is  with  the  models.  In  the  MRAC  system  the  model  represents  the  ideal  plant  response  to 
the  input.  The  model  in  the  MBC  system  is  a  simulation  of  the  plant.  The  adaptation  in 
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the  MBC  system  occurs  in  the  simulation  by  adjusting  the  parameters  of  the  simulation 
through  an  adaptation  law  (tracking  filter)  so  that  the  outputs  of  the  simulation  match 
certain  key  outputs  of  the  plant  when  subjected  to  the  same  inputs.  The  simulation  can 
then  compute  the  desired  feedback  values,  once  the  error  between  the  outputs  of  the  plant 
and  simulation  is  below  some  threshold. 

In  addition,  the  idea  of  robustness  became  paramount  when  working  with  the 
turbine  engine  and  developing  a  control  design  philosophy.  The  engine  behaves  or 
responds  differently  not  only  at  the  many  flight  conditions  at  which  it  will  operate,  but 
also  the  characteristics  will  be  different  based  upon  the  throttle  setting  at  a  given  flight 
condition.  Capturing  of  these  differences  and  including  them  in  a  single  design 
methodology  is  a  very  challenging  task.  In  every  step  of  the  design  process,  both 
structured  and  unstructured  uncertainty  is  considered  to  ensure  robustness  of  the  design. 

The  introduction  of  Quantitative  Feedback  Theory  (QFT)  as  a  design  tool  in  the 
model-based  control  structure  is  a  significant  benefit  to  the  control  system  designer.  QFT 
allows  the  establishment  of  constraints  in  an  orderly  and  concise  manner  so  that  a  low- 
order,  robust  controller  is  found  that  satisfies  overall  system  specifications.  Also,  the 
development  presented  in  this  work  establishes  a  solid  mathematical  foundation  for  the 
design  of  robust  model-based  control  system. 

Another  key  area  of  contribution,  which  is  covered  in  Chapter  3,  is  the 
development  of  a  robust  tracking  filter.  The  tracking  filter  is  the  heart  of  any  model- 
based  control  system.  If  the  model  can  not  be  adjusted  so  that  it  is  an  accurate 
representation  of  the  plant  being  controlled,  then  the  variables  computed  by  the  model,  in 
particular  the  ones  that  will  ultimately  be  used  for  feedback,  will  not  be  accurate  enough. 


77 


In  this  work,  a  method  for  updating  the  states  of  the  model  using  proportional-plus- 
integral  observer  theory  is  presented.  This  concept  was  then  extended  to  the 
multivariable  case. 

The  next  area  of  contribution,  presented  in  Chapter  4,  was  a  further  refinement  of 
the  use  of  the  Perron  root  of  the  interaction  matrix  as  a  measure  of  the  generalized 
diagonal  dominance  of  uncertain  multivariable  plants.  This  capitalizes  on  non-negative 
matrix  theory  and  the  use  of  H-  and  M-matrices.  Earlier  works  present  a  method  for 
designing  a  feedback  loop  around  the  plant  in  order  to  achieve  a  certain  level  of 
decoupling  (Nwokah,  et.al.,  1995).  However,  in  this  work  a  methodology  for  designing  a 
dynamic  forward  loop  precompensator  is  developed.  A  static  precompensator  was 
initially  considered  for  this  design.  However,  the  interaction  index  of  the  plant  was  very 
sensitive  to  the  design  of  a  static  precompensator.  This  was  exhibited  by  a  push-up  in  the 
interaction  index  at  low  frequencies,  which  is  undesirable.  By  designing  a  dynamic 
compensator  the  engineer  now  has  the  ability  to  shape  the  interaction  index  curve, 
thereby  allowing  the  interaction  index  to  stay  small  at  low  frequencies.  The  design  of  the 
precompensator  is  not  a  standard  part  of  the  QFT  design  process.  But,  for  any  practical 
controller  design  the  engineer  desires  to  have  a  system  that  is  as  decoupled  as  possible. 
The  inclusion  of  this  block  permits  the  designer  to  attempt  to  achieve  this  desirable 
characteristic.  Then,  using  the  Quantitative  Feedback  Theory  the  overall  control  design 
process  is  greatly  simplified. 

Finally,  in  Chapter  5,  the  above  theoretical  and  practical  developments  were 
applied  to  an  engine  control  problem.  This  problem  is  unique  in  that  the  feedback 
variables  are  not  directly  sensed.  They  are  computed  within  the  embedded  model  and 
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then  used  to  compute  error  signals  for  the  controller.  In  order  to  include  additional 
realism  into  the  problem,  transport  delay  and  combustion  delays  were  considered  in  the 
model  of  the  engine.  Time  delay  was  also  considered  for  the  execution  of  the  control  law 
and  for  the  sensors  and  actuators.  The  inclusion  of  this  additional  time  delay  made  the 
design  work  much  more  difficult  to  accomplish.  However,  the  controller  demonstrated 
exceptional  time  and  frequency  domain  behaviors,  while  requiring  minimal  gain  over  the 
control  bandwidth. 

6.2  Conclusions 

A  number  of  important  lessons  were  learned  during  the  execution  of  the  design. 
The  most  significant  of  these  concerned  the  limitation  of  achievable  system  performance 
in  the  presence  of  system  time  delay.  This  time  delay  caused  many  difficulties  while 
trying  to  shape  the  loops  to  achieve  overall  system  level  specifications.  Methods  for 
dealing  with  non-minimum  phase  systems  are  known.  However,  it  appears  that 
accounting  for  this  behavior  and  its  effects  becomes  a  unique  aspect  of  every  design. 

Another  key  lesson  was  the  push-pop  effect  that  occurs  when  designing  static 
precompensators.  Even  though  this  phenomenon  was  noticed  in  the  design  structure 
presented  by  Nordgren  (1994),  it  appears  to  be  more  pronounced  problem  when  the 
decoupling  system  is  in  the  forward  path.  This  effect  made  it  necessary  to  design  a 
dynamic  compensator  to  keep  the  interaction  between  the  control  loops  reasonable.  This 
design  was  relatively  simple  and  offered  a  lot  of  benefit  to  the  overall  process. 

To  the  best  of  our  knowledge,  there  exists  no  robust,  model-based,  multivariable 
control  design  methodology  for  the  engine  control  problem.  Further  use  of  this  concept 
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and  refinement  of  the  methodologies  will  greatly  improve  the  overall  engine  system. 
This  could  especially  be  true  in  other  areas  related  to  control  such  as  engine  health 
monitoring. 

6.3  Recommendations 

There  are  many  additional  topics  or  areas  of  improvement  that  became  apparent 
during  the  performance  of  the  work  presented  herein.  The  first  concerns  the  modeling  of 
the  engine.  During  an  investigation  of  existing  engine  modeling  techniques,  different 
approaches  for  modeling  the  dynamics  of  the  system  were  identified.  Some  of  the 
models  were  iterative  and  would  converge  to  a  steady  state  answer  before  returning 
values.  While  others  accepted  small  errors  during  the  transient  and  would  provide  values 
through  each  pass  of  the  code.  While  each  process  will  work,  careful  consideration  must 
be  given  to  the  time  delay  and  overall  impact  of  carrying  small  errors  through  a  transient 
will  have  on  the  model-based  control  concept.  Better  methods  of  representing  the 
components  may  offer  a  significant  advantage  of  current  representations.  Typical  engine 
models  use  lookup  tables  to  represent  the  components;  this  is  the  method  used  for  the 
nonlinear  engine  model  in  this  work.  However,  to  improve  the  functionality  of  the  model 
a  neural  net  could  be  used  to  represent  the  components  thereby  possibly  capturing 
additional  component  characteristics. 

The  second  concerns  the  expansion  of  the  design  to  include  the  remainder  of  the 
flight  envelope.  This  will  not  be  a  major  concern  as  long  as  plant  response  characteristics 
are  similar  and  the  plant  models  are  path  connected.  However,  if  this  is  not  the  case  then 
multiple  controller  designs  are  required  to  cover  the  entire  envelope.  If  multiple  designs 


are  required  then  a  methodology  for  connecting  or  translating  between  the  different 
controller  designs  must  be  used.  One  approach  presented  by  Polly,  et.  al.,  (1988),  uses  a 
curve  fitting  approach  to  connect  all  the  gains  of  the  controllers  together.  It  is  not  known 
if  this  approach  will  work  well  with  this  design  technique  or  if  another  approach  must  be 
found. 

While  the  tracking  filter  appears  to  be  such  a  simple  element  in  the  overall 
concept,  the  importance  of  this  element  can  not  be  over  emphasized.  To  date,  only 
tracking  filters  that  work  well  for  a  single  mode  at  a  time  have  been  found.  The  desire  is 
to  find  a  multipurpose  tracking  filter,  a  filter  that  works  well  for  all  modes  of  operation. 
The  basic  limitation  is  the  number  of  sensed  values  from  the  plant  that  can  be  used  to 
update  parameters  in  the  model.  Research  needs  to  continue  in  order  to  find  a  tracking 
filter  that  meets  a  multitude  of  operational  needs  or  identifies  a  satisfactory  approach  to 
having  multiple  tracking  filters. 

The  final  area  of  interest  concerns  mode  selection.  In  this  work  a  single  mode 
(performance)  was  shown.  One  advantage  of  model-based  control  is  that  the  engine  can 
have  more  than  one  mode  of  operation.  Consider  a  fighter  aircraft.  During  the  lifetime 
of  the  aircraft  it  will  be  used  in  many  different  roles  such  as  combat,  training  and  escort 
missions.  Each  of  these  roles  uses  and  requires  different  capabilities  from  the  same 
engine  hardware.  For  combat,  performance  is  of  primary  concern  while  for  training,  the 
protection  and  safety  of  the  hardware  is  more  important  and  for  the  escort  mission,  fuel 
bum  is  of  interest.  Also,  the  pilot  could  switch  between  the  different  modes,  depending 
upon  the  current  mission  need.  The  use  of  the  model-based  control  concept  will  bring 
new  and  powerful  capabilities  to  the  engine  system  designer  and  user  in  the  future. 
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APPENDIX  A 


EQUIVALENCE  RELATION  DERIVATION 

Consider  a  linear  continuous  time-invariant  system  described  by  the  state-space 
equations 


x  =  Ax  +  Bu 

[A.l] 

y  =  Cx  +  Du 

[A.2] 

where  x  is  an  nxl  vector  of  state-variables,  u  is  an  mxl  vector  of  input  functions,  and  y  is 
an  lxl  vector  of  outputs. 

Solving  equation  [A.l]  with  the  initial  condition  vector  x(0),  and  substituting  into 
equation  [A.2],  we  obtain 

t 

y(t)  =  CeA,x(0)  +  J  CeA(t_T)Bu(x)dx  +  Du(x)  [A.3] 

o 

From  [A.3],  x(0)  can  be  determined  given  the  vectors  y(t)  and  u(t)  over  an  interval  if,  and 
only  if,  the  columns  of  the  matrix  CeAt  are  linearly  independent.  Hence,  a  necessary  and 
sufficient  condition  for  the  system  to  be  observable  is  that  the  rank  of  <j>0  =  n,  where  <[>0 
is  the  matrix 


ct,atcMat)V,..,(at)V 


[A.4] 
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A.1  Observability 

The  system  described  by  [A.l]  and  [A.2]  is  said  to  be  observable  if  there  exists  a 
time  ti>0  such  that  given  the  vectors  u  and  y  over  the  interval  (0,ti)  it  is  possible  to 
deduce  the  initial  state-vector  x(0). 

If  <j)0  has  rank  n  then  the  system  is  completely  observable.  If  <j>0  does  not  have 
rank  n,  then  some  of  the  system  states  are  unobservable  and  have  no  influence  on  the 
system  outputs.  The  rank  deficiency  of  <j)0  indicates  the  number  of  unobservable  states 
but  does  not  identify  these. 

A.2  Observable  Canonical  Forms 

If  a  system  [A,  B,  C]  is  observable,  then  the  matrix  <j)0  defined  in  equation  [A.4] 
has  rank  n.  A  nonsingular  nxn  transformation  matrix  T  can  then  be  formed  from  the 
columns  of  <|)0  as  follows: 

(a)  Inspect  the  columns  of  <j>0  from  the  left  to  the  right  and  retain  only  those  vectors 
which  are  linearly  independent  of  those  previously  selected. 

(bj  Arrange  the  n  linearly  independent  vectors  so  selected  to  form  a  new  matrix  (T), 
where 


[A.5] 


The  integers  (p pf)  are  known  as  the  observability  indices  of  the  system  ,  and 

i>=n- 

i=l 
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The  parameter 


u0  =  max  Pi;  i  =  \,l 


[A.6] 


is  often  called  the  observability  index  for  the  system, 
(c)  Let  T'1  be  described  in  terms  of  its  rows  as 


r_1  = 


Yn 


[A.  7] 


and  let  y  [.  be  the  k*  row  of  T  1 ,  where 


ki  =Ziij;  i=1>^ 

j=! 


[A-8] 


(d)  Using  the  vectors  y  [  the  required  transformation  matrix  T  is  formed  as 


T  =  [yki,Ayki,....A'“-,yki,ykl,Ayki,...,A,l‘-1ykt] 


[A.9] 


(e)  The  observable  standard  form  for  the  system  [A,  B,  C]  is  then  given  by 


A  =  T'1  AT 
B  =  T-1B 
C  =  CT 


[A.  10] 


where 
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[A.  11] 


and  the  diagonal  blocks  Au  are  in  the  companion  form 


0  0 
1  0 


0  X 
0  X 
0  X 


0  0  IX 


[A- 12] 


for  i  =  \,t ;  with  dimensions  x  [it .  The  X’s  stand  for  possible  nonzero  entries.  The 
off-diagonal  blocks  Ay,  i  *j,  have  the  form 


X 

X 


X 

X_ 


[A- 13] 


The  transformed  output  matrix  C  is 

C  =  [C1,C2,...,Cf]  [A.  14] 

where  the  l  x  p.;  blocks  C;  have  the  special  form 

Ci=[0  ;  £i]  [A- 15] 

with 
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0  for  j  <  i 
1  for  j  =  i 
X  forj  >  i 


[A.  16] 


s 


i.e.  the  entries  below  the  unit  entry  in  the  last  column  of  C;  are  not,  in  general,  equal  to 
zero.  These  nonzero  entries  can  be  readily  removed,  if  so  desired,  by  a  simple 
nonsingular  transformation  of  the  system  outputs. 

The  matrix  B  has  no  special  form. 


92 


APPENDIX  B 


NECESSARY  PROOFS 


B.l  Proof  of  Theorem  4.1 

Without  loss  of  generality,  assume  that  Xz  <  0.5.  We  then  need  to  show  that 
Xt<\  and  hence  Z  <=H.  Consider  the  diagonal  regular  splitting  of  Z  (this  is  always 
possible  whenever  Z  e  H  )  as: 

Z  =  D  +  C  =  (l  +  CD-1  )D  =  (I  +  M)D  e H  [B.l] 

Note  that  I  +  M  6  H  since  multiplication  by  a  nonsingular  complex  diagonal  matrix  does 
not  destroy  the  H-matrix  characteristic.  Since 

p(m)sp(m.)=p(c.d;,)->.z<i  [b.2] 


where  p(»)  denotes  the  spectral  radius  of  •,  we  may  write  Z  as  the  following  convergent 
series: 


Z  =  D-1  (I  +  M)-1  =  D_1 


f  oo  ^ 

k=l  J 


[B.3] 


As  D  is  diagonal,  ZeHoI  +  ^](-M)k  €  H ,  or  equivalently,  if  and  only  if: 

k=l 
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p(i(-M)k)  =^<i- 


[B.4] 


Now, 


f  “  .  \ 


^=P 


co  \  oo 


2(-M)1  spix  hZpK)=S^=rr-  [b-5] 

Vk=l  J+  \k=l  J+  k=l  k=l  1 


Then  the  infinite  series  sum  [B.5]  reduces  to: 


„  Xz  X 
X^< — —  o 


z  -  "z  <x.. 


z  \-xz  \+xt 


[B.6] 


Thus, 


A,z  <  0.5  =>  <  1  =?>  Z  e  H 


[B.7] 


which  completes  the  proof. 

B.2  Proof  of  Theorem  4.3 

Given  that  the  elements  K  (precompensator),  G  (diagonal  controller)  and  F 
(prefilter)  are  stable.  Then  TY/R  is  unstable  if  and  only  if 

det(Z)  =  det((l  +  M)(Pd  +  G))  =  det(I  +  M)det(PD  +  G)  [B.8] 


has  zeros  in  the  right-hand  plane. 
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Thedet(l  +  M)has  no  zeros  in  the  right  hand  plane  if  Xp(M)  <  1.  In  addition, 


the  det(PD  +  G)  = 

i=l 


has  no  zeros  in  the  right  hand  plane  if  each  (l  +  g^jj) 


has  no  zeros  in  the  right-hand  plane.  This  completes  the  proof. 
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APPENDIX  C 


COMPUTER  ROUTINES 


Contained  herein  is  a  set  of  routines  written  for  Matlab  4.0  that  were  used  to 
generate  the  design  given  in  Chapter  5.  These  routines  can  be  grouped  into  categories  for 
design,  analysis  and  support.  Below  is  a  listing  of  each  routine  by  category  and  the  page 
on  which  it  may  be  found. 


Design  Page  # 


a)  decouple .  97 

b)  gdesl. .  99 

c)  gdes2 . 101 

d)  pfdesl . 103 

e)  pfdes2 . 105 


Analysis 


a)  interact . 107 

b)  mybodes . 109 

c)  my  steps . ■ . Ill 

d)  pathconn . . . . . 114 

Support 

a)  abed . 116 

b)  abcd_n . 117 

c)  actblck . 118 

d)  actblck_w . 120 

e)  gk . 122 

f)  gk_w . 124 

g)  pfilter . 126 

h)  pfilter  w . 127 

i)  perron . 128 
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% 

%  Name  :  Decouple 


%  Purpose:  To  generate  data  for  the  design  of  the 
%  precompensator 


%  Inputs  :  None  passed  to  the  routine  must  have  system 
%  matrices  available  and  an  appropriate  frequency 

%  vector  must  be  set  within  the  routine 

% 

%  Outputs:  (qij/qjj)  needed  for  loop  shaping 
%  bounds  for  loopshaping 

% 

%  Loopshaping  Commands  : 

%  element  12  :  lpshape (w,bndl2,-ql2_data (1,  : ) ,  []  ) 

%  element  21  :  Ipshape (w, bnd21, -q21_data (1, : ) , [ ] ) 


%  Written  by  :  Zane  D.  Gastineau 
% 

%=================“““==“== 


%  Load  system  matrices 

load  c: \matlab\turbojet\normmat\sls_an 
load  c:\matlab\turbojet\normmat\sls_bn 
load  c: \matlab\turbojet\normmat\sls_cn 
load  c:\matlab\turbojet\normmat\sls_dn 

%  Frequency  range  of  concern 

w=sort ( [logspace (-1,2, 50), 0.5, 1.0, 1.5, 3. 0,6. 0,8.0]); 


%  Housekeeping 
%=======t===========; 

j=sqrt (-1) ; 
q21_data= []  ; 
ql2_data=[]  ; 
v=0; 


%  ====================—========="=:= ====—— ~===========:::::— — 

%  Build  up  plant  and  generate  loop  shaping  data. 
%—~~===——~——~~====—~—~===———=====———=====z——:::==::=:—=:=:=:==~=~ 

for  k=15 : length (an) -2 
k 

v=v+ 1 ; 

% 

[ap, bp, cp, dp] =abcd_n (an, bn, cn, dn,  k)  ; 

% 

for  mi=l : length (w) 

P=cp*  (inv  ( j  *w  (mi)  *eye  (size  (ap)  )  -ap)  *bp)  +dp; 
[act, senl^actblck^wfwtmi) ) ; 

L=sen*P*act ; 
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Lhat=inv (L) ; 

ql2_data (v,mi) = ( 1/Lhat (1, 2) ) *Lhat (2,  2) ; 
q21_data (v,mi)=(l/Lhat (2, 1) ) *Lhat (1, 1) ; 

end 

end 


%  Generate  Bounds  for  the  precompensator 
% 

%  | 1-kij (qij/qj j ) I  <=  1 

%  Frequencies  to  Generate  bounds  at 
% 

w_bnd= [0.1  0.51368]; 

bndl2=genbnds  (10,  w,  w_bnd,  1, 1,  -ql2_data,  1,  0,  -ql2__data  (1,  : )  ) 
bnd21=genbnds (10, w, w_bnd, 1, 1, -q21_data,  1, 0, -q21_data (1, : ) ) 
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o\P  o\°  o\o  o\°  o\o  o\o  o\o  o\P  o\°  0\o  o\P  o\o  o\P  o\o  <A° 


% 

%  Name  :  gdesl.m 
% 

Purpose  :  Routine  to  compute  loop  1  bounds 

Inputs  :  None  passed  to  the  routine  must  have  system 

matrices  available  and  an  appropriate  frequency 
vector  must  be  set  within  the  routine. 

Outputs  :  Bounds  for  use  in  the  loop  shaping  routine  of  the 
QFT  toolbox. 

Loopshaping  Command  : 

Ipshape (w, bnd, qsll (1, : ) ,  [] ) 

Written  by  :  Zane  D.  Gastineau 


%  Load  plant  set 

load  c :  \mat lab\turbo j  et  \normmat  \sls__an 
load  c:  \mat lab\ turbo j  et  \normmat\sls__bn 
load  c: \matlab\turbojet\normmat\sls_cn 
load  c: \matlab\turbojet\normmat\sls_dn 


%  Load  Frequency  set  points 

w=sort { [logspace (-1,2,50) ,0.8, 1.0, 3.0, 6.0,8.0,10,12,20,4.4984] ) ; 


%  Housekeeping 

i=sqrt (-1) ; 
v=0; 

%  Build  up  Plant  over  the  frequency  range 

for  p=15 : length (an) -2 
v=v+l; 

[ap,bp,cp,dp]=abcd_n(an,bn,  cn,dn,p)  ; 
for  k=l: length (w) 

P=cp*inv (w (k) *i*eye (size (ap) ) -ap) *bp+dp; 

[G, KK] =gk_w ( w ( k) ) ; 

[Act,Sen]=actblck_w(w(k)); 

L=S  en*  P  * Act  *  KK ; 

lsll (v,  k) =L (1, 1) ;lsl2 (v, k) =L (1, 2) ; 
ls21(v,k)=L(2,l) ; ls22 ( v, k) =L ( 2, 2 ) ; 

PKhat=inv (L) ; 

qsll (v, k)  =  1/PKhat ( 1, 1 ) ; qsl2 (v, k)  -  1/PKhat (1, 2) ; 
qs21 (v, k)  =  1/PKhat (2, 1) ;qs22 (v, k)  =  1/PKhat (2, 2 ) ; 

end 
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end 


%  Set  and  generate  Bounds  for  QFT  design 

%  Compute  Upper  (emu)  and  Lower  (cml)  bound  over  frequency 

cmul=freqcp (64, [1  9.6  64],w); 
cmu2=freqcp (9, [1  3.6  9], w) ; 
cmu=max ( emul , cmu2 ) ; 
cmll=freqcp (169, [1  31.2  169],w); 
cml2=freqcp (36, [1  14.4  36], w); 
cml=min (cmll, cml2) ; 

Ws=abs ( [emu;  cml] ) ; 

bndl  =  sisobnds (7, w,  [0 . 8  1.0  3.0  6.0] ,Ws,qsll,  []  ,  2)  ; 
bnd3  =  sisobnds (2, w,  [6] ,  1/0. 4352, qsll,  [] , 2) ; 
bnd4  =  sisobnds (2, w,  [8] ,  1/0 . 51264, qsll,  [], 2) ; 
bnd=grpbnds (bndl,bnd3,bnd4) ;plotbnds (bnd) 
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% 

%  Name  :  gdes2.m 

% 

%  Purpose  :  Routine  to  compute  loop  2  bounds 
% 

%  Inputs  :  None  passed  to  the  routine  must  have  system 
%  matrices  available  and  an  appropriate  frequency 

%  vector  must  be  set  within  the  routine. 

q_ 

'o 

%  Outputs  :  Bounds  for  use  in  the  loop  shaping  routine  of  the 
%  QFT  toolbox. 

o, 

"o 

%  Loopshaping  Command  : 

%  lpshape (w,bnd, qs22 (1,  : )  ,  [] ) 

o. 

"o 

%  Written  by  :  Zane  D.  Gastineau 
% 


ii 
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it  i 

load  c: \matlab\turbojet\normmat\sls_an 
load  c:\matlab\turbojet\normmat\sls_bn 
load  c:\matlab\turbojet\normmat\sls_cn 
load  c:\matlab\turbojet\normmat\sls_dn 

Q _ _ _ — — — ssssaass=====s==sassi=aB=e====SSS===as— : ===== ==: ===== 

"O - 

%  Load  Frequency  'set  points 

w=sort ( [logspace (-1,2,50) ,0.8,1.0,3.0,6.0,8,10,12,20] ) ; 


%=========================================== 

%  Housekeeping 

i=sqrt (-1) ; 
v=0; 

%  Build  up  plant  over  the  frequency  range 


for  p=15 : length (an) -2 
v=v+l 

[ap,  bp,  cp,  dp]  =abcd_n  (an,  bn,  cn,  dn,p)  ; 
for  k=l: length (w) 

P=cp*inv(w(k) *i*eye (size (ap) ) -ap) *bp+dp; 

[G, KK] =gk_w (w(k) ) ; 

(Act,Sen]=actblck_w(w(k) )  ; 

L=Sen*P*Act*KK; 

lsll (v, k) =L (1, 1) ; Is 12 (v, k)=L(l,2) ; 
ls21 (v,k)=L(2,l) ;ls22(v,k)=L(2,2) ; 

PKhat=inv (L) ; 

qsll (v, k)  =  1/PKhat (1, 1) ;qsl2 (v,  k)  =  1/PKhat (1,  2)  ; 
qs21 (v, k)  =  1/PKhat (2, 1) ;qs22 (v, k)  =  1/PKhat (2,  2 )  ; 

end 
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end 


%  Set  and  generate  Bounds  for  QFT  design 

%  Compute  Upper  (emu)  and  Lower  (cml)  bound  over  frequency 
%  =——~~====——~==— 
cmu=freqcp (36,  [1  60  36], w); 
cml=freqcp (16, [1  40  16], w) ; 

Ws=abs ( [emu;  cml] ) ; 

bndl  —  sisobnds (7 , w, [0 . 8  1.0  3.0  6. 0] ,Ws, qs22, [] , 2) ; 
bnd3  =  sisobnds (2, w,  [6] ,  1/0 . 4352, qs22,  [ ] , 2 ) ; 
bnd4  =  sisobnds (2, w, [8], 1/0. 5 12 6, qs 22, [ ] , 2) ; 
bnd=grpbnds {bndl , bnd3 , bnd4 ) ; 
plotbnds (bnd) 
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%  Name  :  pfdesl.m 

% 

%  Purpose  :  To  generate  required  data  for  the  design  of  the 
%  prefilter 

% 

%  Inputs  :  None  passed  to  the  routine  must  have  system  matrices 
%  available  and  an  appropriate  frequency  vector  must 

%  be  set  within  the  routine. 

% 

%  Outputs  :  L(i,j)  for  prefilter  design  loopshaping. 

% 


%  Loopshaping  Command  : 

%  pf shape (7, w, w, Ws, Lll  (1, : ) , [] , Gil (1, : ) ) 

% 

%  Written  by  :  Zane  D.  Gastineau 


load  c: \matlab\turbojet\normmat\sls_an 
load  c: \matlab\turbojet\normmat\sls_bn 
load  c: \matlab\turbojet\normmat\sls_cn 
load  c: \matlab\turbojet\normmat\sls_dn 


Load  Frequency  set  points 


w=sort ( [logspace (-1, 2, 50), 0.8, 1.0, 3. 0,6. 0,10,  12,20,4.4  984]); 


%  Housekeeping 


Build  up  Plant  over  the  frequency  range 


for  p=15 : length (an) -2 
v=v+l 

[ap,bp,  cp,  dp]  -abcd__n  (an,  bn,  cn,  dn,p)  ; 
for  k=l : length (w) 

P=cp*inv (w (k) *i*eye (size (ap) ) -ap) *bp+dp; 
[G, KK] =gk_w ( w ( k ) )  ; 

Gil (v, k) =G (1, 1) ; 

[Act, Sen] =actblck_w (w ( k) ) ; 
L=Sen*P*Act*KK; 

Lll ( v, k) =L ( 1, 1 ) ; L12 ( v, k) =L ( 1, 2) ; 

L21  (v,  k)  -L  (2, 1)  ; L22  (v,  k)  =L  (2, 2 )  r 

end 

end 
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%  Compute  Upper  (emu)  and  Lower  (cml)  bound  over  frequency 


cmul=freqcp(64, [1  9.6  64],w); 
cmu2=freqcp ( 9, [1  3.6  9],w); 
cmu=max  ( emul ,  cmu2 )  ; 
cmll=freqcp (169, [1  31.2  169], w); 
cml2=freqcp (36, [1  14.4  36], w); 
cml=min ( cml 1 , cml 2 ) ; 

Ws=abs ( [emu;  cml] ) ; 
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%  Name  :  pfdes2.m 
% 

%  Purpose  :  To  generate  required  data  for  the  design  of  the 
%  prefilter 

% 

%  Inputs  :  None  passed  to  the  routine  must  have  system  matrices 
%  available  and  an  appropriate  frequency  vector  must 

%  be  set  within  the  routine. 

% 

%  Outputs  :  L (i, j )  for  prefilter  design  loopshaping. 

% 

%  Loopshaping  Command  : 

%  pf shape (7, w, w, Ws, L22 (1, : ) , [] , G22 (1, : ) ) 

% 

%  Written  by  :  Zane  D.  Gastineau 

% 


H 

ll 

ii 

ll 

ii 

ll 

ll 

il 

1 

i 

Load  plant  set 

load  c:\matlab\turbojet\normmat\sls_an 
load  c:\matlab\turbojet\normmat\sls_bn 
load  c: \matlab\turbojet\normmat\sls_cn 
load  c:\matlab\turbojet\normmat\sls_dn 

%=========================================================== 

%  Load  Frequency  set  points 

w=sort ( [logspace (-1, 2,50) ,0.8, 1.0, 3.0, 6.0, 10, 12,20,4.4984] ) ; 


%============== 

%  Housekeeping 
%============== 

i=sqrt (-1) ; 
v=0; 


%  Build  up  Plant  over  the  frequency  range 

for  p=15 : length (an) -2 
v=v+l 

[ap,  bp,  cp,  dp]  =abcd_n  (an, bn,  cn,  dn,  p) ; 
for  k=l: length (w) 

P=cp*inv(w(k) *i*eye (size (ap) ) -ap) *bp+dp; 
[G, KK] =gk_w (w(k) ) ; 

G22 (v, k) =G (2, 2) ; 

[Act ,  Sen] =actblck_w ( w ( k) )  ; 

L=Sen*P*Act*KK; 

Lll (v, k)=L(l, 1) ;L12 (v,k)=L(l, 2) ; 

L21 (v, k) =L (2, 1) ;L22 (v,k)=L(2,2) ; 

end 

end 
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%  Compute  Upper  (emu)  and  Lower  (cml)  bound  over  frequency 


cmu=freqcp (16, [1  16  16] ,  w); 
cml=freqcp (36, [1  24  36], w) ; 
Ws=abs ( [emu;  cml] ) ; 
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%=_^=======  ================:======  =  ^==:==========:====:============: 

% 

%  Name  :  Interact. m 
% 

%  Purpose  :  Routine  to  compute  the  interaction  index 
%  for  the  system 

% 

%  Inputs  :  None  passed  to  the  routine  must  have  system 
%  matrices  available  and  an  appropriate  frequency 

%  vector  must  be  set  within  the  routine, 

% 

%  Outputs  :  Plot  of  the  interaction  index  for  either  the 
%  direct  or  for  the  inverse  system. 

o. 

~o 

%  Written  by  :  Zane  D.  Gastineau 
% 

%=================================  =  =  -===-==:================= 


%  Load  plant  set 

load  c:\matlab\turbojet\norimat\slsjan 
load  c: \matlab\turboj et \normmat \sls_bn 
load  c:\matlab\turbojet\normmat\sls_cn 
load  c: \matlab\turboj et\normmat\sls_dn 


%  Housekeeping 
%=====  =========: 

v=0; 

i-sqrt (-1) ; 


%  ======================= 

%  Load  Frequency  Data 


w=sort ( [logspace {-2, 3, 50), 0.5, 1.0,3. 0,6. 0,8.0] ) ; 


% - 

%  Compute  interaction  index 


for  k=15 : size (an, 1) -2 
k 


v=v+l; 

[ac, be, cc, dc] =abcd_n (an, bn,  cn,  dn,  k) ; 
for  q=l : length (w) 

[act, sen] =actblck_w (w (q) ) ; 

%ke=eye (2) ;  %  Use  this  for  plant  interaction  index 

[G, ke]=gk_w(w(q) ) ; 

P=cc*inv (w (q) *i*eye (3) -ac) *bc+dc; 
vl=sen*P*act*ke; 

%vl=inv (P) ;  %Use  this  for  the  inverse  plant. 
vld=diag (diag (vl) ) ; 

vlo-vl-vld; 

Ml=vlo*inv (vld)  ; 

[lambda, v0]=perron (Ml) ; 
zl (v, q) =lambda; 
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end 


end 

%  Plot  Interaction  Index 
for  q=l:v 

semilogx (w, zl (q,  : )  ) 
hold  on 

end 

grid 

xlabel ( 1  frequency  (rad/s) 1 ) 
ylabel ( f Interaction  Index1) 
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%  Name  :  mybodes .m 

% 

%  Pupose  :  To  plot  magnitude  vs  frequency  for  different 
%  functions. 

% 

%  Inputs  :  None  passed  to  the  routine  must  have  system 
%  matrices  available  and  an  appropriate  frequency 

%  vector  set  within  the  routine.  Actuator, 

%  Sensor,  Controller  and  Precompensator  must  be 

%  defined  outside  the  routine  to  be  called 

% 

%  Outputs  :  Magnitude  responses  for  the  Plant,  Plant + 

%  precompensator,  Plant tprecompensator+cont roller, 

%  closed  loop  system  and  closed  loop  system  with 

%  prefilter. 

% 

%  Written  by  :  Zane  D.  Gastineau 
% 


%  Select  which  plot  you  want. 


dispC Select  the  desired  plot:1) 
disp  ( '*  1 ) 
dispC  1. 
disp{*  2. 
disp(T  3. 
disp(?  4. 
disp{*  5. 
disp ( 1  ') 


in  -  input ( 1 


Plant  ( P) ' ) 

Plant -HPrecompensator  (PK)  1 ) 
Plant+Precompensator+Controller  (L) 1 ) 

Closed  loop  system  [inv (I+L) *L] f ) 

Closed  loop  system  with  prefilter  [inv (I+L) *LF] 1 ) 
Enter  choice  :  1 ) ; 


%  Load  plant  set 

%======—= ========-======-==-=========- 

load  c:  \matlab\turbojet\normmat\sls_an 
load  c:\matlab\turbojet\normmat\sls_bn 
load  c:  \matlab\turbojet\normmat\sls_cn 
load  c:  \matlab\turboj  et\normmat \sls__dn 

_ _ _ 

"O - — - 

%  Load  Frequency  set  points 


w=sort ( [logspace (-1,2,50) ,0.8,1.0,3.0,6.0,10,12,20,4.4984] ) ; 
%  Build  up  Plant  over  the  frequency  range 

Cl _ _ _ — - = 

"O - - - 

for  p=15 : length (an) -2 
P 

[ap,bp,cp,dp]=abcd_n(an,bn,cn,dn,p)  ; 
for  k=l: length (w) 

Plant=cp*inv(w(k) *i*eye (size (ap) )-ap) *bp+dp; 


o\o  <*o  c\o  o\o 


[Act, Sen] =actblck_w (w ( k) ) ; 

P=Sen*Plant*Act ; 

Pll (p,  k)  =P  (1, 1) ; P12 (p, k)=P(l,  2) ; 

P21 (p, k) =P (2, 1) ; P22 (p, k) =P (2, 2) ; 

[G,  KK] =gk_w (w (k) ) ; 

PK=P*KK; 

PK11 (p,  k) =PK(1, 1) ; PK12 (p, k) =PK (1, 2 ) ; 

PK21 (p, k) =PK(2, 1) ; PK22 (p, k)=PK(2,2) ; 
PKhat=inv (PK) ; 

qsll (p, k) =l/PKhat (1,1) ;qsl2 (p,  k) =l/PKhat (1, 2) ; 
qs21 (p, k) =l/PKhat (2,1) ; qs22 (p, k) =l/PKhat (2, 2) ; 
L=PK*G; 

Lll (p, k) =L (1, 1) ;L12(p,k)-L(l,2)  ; 

L21  (p,  k)  =L  (2, 1)  ; L22  (p,  k)  =L  (2 , 2 )  ; 

T=inv (eye (2) +L) *L; 

Til  (p,  k)  =T  (1, 1)  ;  T12  (p,  k )  =T  (1,2)  ; 

T21 (p, k) =T (2, 1) ;T22 (p, k)=T (2, 2) ; 

[F] =pfilter_w(w(k) )  ; 

TF=T*F; 

TF11 (p, k) =TF(1, 1) ;TF12 (p,  k)=TF(l,  2) ; 

TF21 (p, k) =TF (2, 1) ;TF22 (p, k) =TF (2, 2) ; 
end 

end 

Generate  Magnitude  vs  frequency  plots 


compute  individual  magnitude  responses 


for  k=15:29 
if  in  ==  1 

subplot (2,2,1) ; semilogx (w, db (Pll (k, : ) ) ) ; grid; hold  on; 
subplot (2,2,2) ; semilogx (w, db (P12 (k, : ) ) ) ; grid; hold  on; 
subplot (2, 2, 3) ; semilogx ( w, db ( P2 1 (k, : ) ) ) ; grid; hold  on; 
subplot (2, 2, 4) ; semilogx (w, db (P22 (k, : ) ) ) ; grid; hold  on; 
elseif  in  ==  2 

subplot (2,2,1) ; semilogx (w, db { PK11 (k, : ) ) ) ; grid;  hold  on; 
subplot (2,2,2) ; semilogx (w, db (PK12 (k, : ) ) ) ; grid; hold  on; 
subplot (2, 2, 3) ; semilogx (w, db ( PK21 (k, : ) ) ) ;grid;hold  on; 
subplot (2,2,4) ; semilogx (w, db ( PK22 (k, ; ) ) ) ; grid; hold  on; 
elseif  in  ==  3 

subplot (2,2,1) ; semilogx (w,db (Lll (k, : ) ) ) ; grid; hold  on; 
subplot (2, 2,2) ; semilogx (w,db (L12 (k, : ) ) ) ;grid;hold  on; 
subplot (2,2,3) ; semilogx (w, db (L21 (k, : ) ) ) ; grid; hold  on; 
subplot (2, 2,4) ; semilogx (w,db(L22 (k, : ) ) ) ; grid; hold  on; 
elseif  in  ==  4 

subplot (2,2,1) ; semilogx (w,db (Til (k, : ) ) ) ; grid; hold  on; 
subplot (2,2,2) ; semilogx (w,db(T12(k,  :) ) ) ; grid; hold  on; 
subplot (2,2,3) ; semilogx (w,db(T21 (k, : ) ) ) ; grid; hold  on; 
subplot  (2,2,4)';  semilogx  (w,  db  (T22  (k,  : ) ) )  ;  grid;  hold  on; 
elseif  in  ==  5 

subplot (2,2,1) ; semilogx ( w, db ( TF11 (k, : ) ) )  ;  grid; hold  on; 
subplot (2, 2, 2) ; semilogx (w, db ( TF12 (k, :) ) ) ;grid;hold  on; 
subplot (2,2,3) ; semilogx (w, db (TF21 (k, : ) ) ) ; grid; hold  on; 
subplot (2, 2, 4) ; semilogx (w, db (TF22 (k, : ) ) ) ;grid;hold  on; 

end 

end 
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% 


% 

%  Name  :  mysteps.m 
% 

%  Purpose  :  To  generate  time  response  plots  of  the  plant, 

%  the  closed  loop  system  and  the  closed  loop 

%  system  with  the  prefilter. 

% 

%  Inputs  :  None  passed  to  the  routine.  Must  have  the 
%  actuator  and  sensor  routines  predefined. 

% 

%  Outputs  :  Plot  of  plant  and  closed  loop  system  step 
%  responses 

% 

%  Written  by  :  Zane  D.  Gastineau 
% 


%  Load  system  matrices 

load  c: \matlab\turboj et\normmat \sls_an 
load  c: \matlab\turbo j et \normmat\sls_bn 
load  c: \matlab\turbojet\normmat\sls_cn 
load  c:  \matlab\turboj  et\normmat \sls__dn 


%  Select  proper  plot  routine 

disp{' Select  the  desired  Plot:') 
disp  ( 1  1  ■) 

disp ( 1  '  1.  Plant  (P) 1 ) 

disp('  2.  Closed  Loop  System  [inv (I+L) *L] ' ) 

disp(’  3.  Closed  Loop  System  with  Prefilter  [inv (I+L) *L*F] 1 ) 

disp(T  ') 

in  =  input ( 1  Enter  Choice  :  1 ) ; 


%  Housekeeping 

hold  off 
cl  f 

%  Load  appropriate  State  Space  Matrices 

[act, sen]=actblck; 

[G, KK] =gk; 

[F] =pfilter; 

[Aact ,  Bact ,  Cact ,  Dact  ]  =branch  ( act )  ; 

[Asen, Bsen, Csen, Dsen] =branch ( sen) ; 

[Ag, Bg, Cg, Dg] ^branch (G) ; 

[Ak, Bk, Ck, Dk] =branch (KK) ; 

[Fa, Fb, Fc, Fd] =branch (F) ; 


111 


o _ 
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%  Build  up  System  Matrices 

%=======================================—————————— 

for  p=16: size (an, 1) -2 
P 

[A, B, C, D] =abcd_n (an, bn, cn,dn,p); 

[A,  B,  C,  D]  =series  (Aact,  Bact,  Cact,  Dact,  A,  B,  C,  D) ; 

[A, B, C, D] =series (A, B, C, D, Asen, Bsen, Csen,  Dsen)  ; 
if  in  ==  1 
for  m=l : 2 
t— [0:0. 1:5]; 

[y] -step  (A,  B,  C,  D,m,  t )  ; 
for  n=l:2 

str- ['22'  num2str (m+ (n-1) *2) ] ; 

.subplot (str) ;grid; 
plot (t, y ( : , n) ) 
hold  on 

end 

end 

subplot (221) ; title ( 'Delta  Fuel  Flow' ); ylabel (' Delta 

Thrust ' ) ; 

subplot (222) /title ( 'Delta  Nozzle  Area'); 

subplot (223) ; ylabel ( 'Delta  Stall  Margin' ) ;xlabel ( 'time 

(sec) '); 

subplot (224) ;xlabel ( 'time  (sec)'); 
elseif  in  ==  2 

[A, B, C, D] =series (Ak, Bk,  Ck,  Dk, A,  B,  C,  D)  ; 

[A, B, C, D] =series (Ag, Bg, Cg, Dg, A, B, C, D) ; 

[A,  B,  C,  D]  =cloop  (A,  B,  C,  D,  -1 )  ; 
for  m=l : 2 

t= [0:0.1:10] ; 

[y] =step (A, B, C, D,m, t ) ; 
for  n-1: 2 

str=['22'  num2str (m+ (n-1) *2) ]  ; 
subplot (str) ; grid; 
plot  (t,  y  (•:  ,n)  ) 
hold  on 

end 

end 

subplot  (221) ;title ( 'Delta  Thrust  Demand'); 

ylabel ( ' Delta  Thrust ' ) ; 

subplot (222) ; title ( 'Delta  Stall  Margin  Demand'); 
subplot (223) ;ylabel ( 'Delta  Stall  Margin'); 

xlabeM'time  (sec)'); 
subplot (224) ;xlabel ( 'time  (sec) ' ) ; 
elseif  in  ==  3 

[A,  B,  C, D] =series (Ak, Bk, Ck, Dk, A, B, C,  D) ; 

[A, B, C, D] -series ( Ag , Bg , Cg , Dg ,A,B,C,D); 

[A,  B,  C,  D]  =cloop  (A,  B,  C,  D,  -1)  ; 

[A,  B,  C,  D]  -series  (Fa,  Fb,  Fc,  Fd,  A,  B,  C,  D)  ; 
for  m=l : 2 

t— [0:0.1:10] ; 

[y] -step (A, B, C, D,m, t) ; 
for  n— 1:2 

str- ['22'  num2str (m+ (n-1) *2) ] ; 
subplot (str) ;grid; 


112 


end 


plot (t, y ( : , n)  ) 
hold  on 


end 

subplot (221) ; title (' Delta  Thrust  Demand'); 

ylabel ( ' Delta  Thrust ’ ) ; 

subplot (222) /title ( 'Delta  Stall  Margin  Demand') 
subplot (223) /ylabel ( 'Delta  Stall  Margin')/ 
xlabel('time  (sec)')/ 
subplot (224) /xlabel  ( 'time  (sec)')/ 

end 

end 
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% 

%  Name  :  Pathconn.m 
% 

%  Purpose  :  Routine  for  investigating  the  path  connectedness 
%  of  one  plant  to  another. 

% 

%  Inputs  :  None  supplied  to  the  routine.  The  user  will  be 
%  prompted  for  the  nominal  plant  # .  Then  the  user 

%  can  test  any  other  plant  within  the  plant  set 

%  one  at  a  time. 

% 

%  Outputs  :  Plot  of  ratio  of  determinants.  Indicates  path 
%  connectedness. 

% 

%  Written  by  :  Zane  D.  Gastineau 
%  6909  Custer  Road  #2104 

%  Plano,  Tx  75023 

%  (972)491-3776 

%  load  state-space  models 

load  c: \matlab\turboj et\normmat \sls_an 
load  c:\matlab\turbojet\normmat\sls_bn 
load  c: \matlab\turbojet\normmat\sls_cn 
load  c:  \matlab\turbojet\normmat \sls__dn 


%=============================:=========:=============:=== 

%  Load  the  nominal  plant 

ini  =  input ('Enter  the  number  for  the  nominal  plant  :  '); 

[anom, bnom, cnom, dnom] =abcd_n (an, bn, cn, dn,  ini) ; 


%  Load  desired  frequency  range 
%===="““==“====“=====s=====s=! 
w=logspace (-1,2,50); 


%  Select  display  method 
disp ( '  ' ) 

disp(' Select  Method  of  Display1) 
disp('  ') 

disp(!  1.  User  defined  plant') 
dispC  2.  Loop  through  all  plants') 
disp(*  ') 

in2  =  input ( '  Enter  Method  :  ' ) ; 

if  in2  ==  1 
disp  ( '  ' ) 

in3  =  input ('Enter  Desired  Plant  to  Check  :  '); 

z=[]  ; 

[ac, be, cc, dc] =abcd_n (an, bn, cn, dn, k) ; 
for  q=l : length (w) 
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vl=cnom*inv  (w  (q)  *i*eye  (3)  -anom)  *bnom+dnom; 
v2=cc*in v (w (q) *i*eye (3) -ac) *bc+dc; 
rl=det (vl) ; 
r2=det (v2) ; 
z=[z;r2/rl] ; 

end 

for  q=l : length (w) 

plot (real (z) , imag (z)  ) 
str=[’ plant  1  num2str (k)  ]  ; 
title (str) 
ylabel ( 1 imag(z) 1 ) 
xlabel ( 1  real ( z)  ’ ) 

end 

elseif  in2  ==  2 

for  k=l : size (an, 1)  %use  to  loop  through  plants 

k 

z«[]; 

[ac,bc, cc, dc] =abcd_n (an, bn, cn, dn,  k) ; 

for  q-1 : length (w) 

vl=cnom*inv (w (q) *i*eye (3) -anom) *bnom+dnom; 
v2=cc*inv (w (q) *i*eye (3) -ac) *bc+dc; 
rl=det (vl) ; 
r2=det (v2) ; 
z=[z;r2/rl] ; 

end 

for  q=l: length (w) 

plot (real (z) , imag(z)  ) 
str=[ ’plant  ’  num2str (k) ] ; 
title (str) 
ylabel ( ’ imag ( z ) ’ ) 
xlabel ( ’ real (z) ’ ) 

end 


end 


pause 


end 
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0\o  o\o  o\o  o\o  o\o  o\o  oY3  o\°  oY>  o\P  d\°  o\P  o\o  oY>  o\°  o\° 


function  [A,  B,  C,  D]  =abcd_n  (a,  b,  c,  d,  k) 


Name  :  abcd.m 

Purpose  :  Function  returns  kth  plant  of  full  set  of  plants. 
All  outputs  available . 

Inputs  :  Must  supply  system  matrices  as  large  matrices  where 
each  plant  makes  up  one  row.  User  must  also  supply 
which  plant  they  are  interested  in. 

Outputs  :  The  desired  plant  in  State  Space  Form 

Written  by  :  Zane  D.  Gastineau 


A  =  [a (k, 1)  a (k, 2)  a(k,3) 
a(k,4)  a(k,5)  a(k,6) 
a (k, 7)  a (k, 8)  a(k,9)]; 

B  =  [  b (k,  1)  b(k,2)  b(k,3)  b(k,4)  b(k,5) 

b  (k,  6)  b  (k,  7)  b(k,8)  b(k,9)  b(k,10) 

b (k, 11)  b (k, 12)  b (k, 13)  b(k,14)  b(k,15)]; 

C  =  [  c  (k,  1)  c  (k,  2)  c  (k,  3) 
c (k, 4 )  c (k, 5)  c (k,  6) 
c (k, 7)  c(k, 8)  c  (k, 9) 
c (k, 10)  c (k, 11)  c (k, 12) 

c (k, 13)  c (k, 14)  c (k, 15) 

c (k, 16)  c (k, 17)  c (k, 18) 

c (k, 19)  c (k, 20)  c (k, 21) 

c (k, 22)  c (k, 23)  c(k,24) 

c (k, 25)  c (k, 26)  c(k,27)]; 

D  =  [  d(k, 1)  d(k,2)  d(k,3)  d(k,4)  d(k,5) 

d(k,  6)  d (k,  7)  d(k,8)  d(k,9)  d(k,10) 

d(k, 11)  d(k, 12)  d (k, 13)  d(k,14)  d(k,15) 

d(k, 16)  d(k, 17)  d(k,18)  d(k,19)  d(k,20) 

d(k, 21)  d(k, 22)  d(k,23)  d(k,24)  d(k,25) 

d(k, 26)  d (k, 27 )  d(k,28)  d(k,29)  d(k,30) 

d(k, 31)  d(k, 32)  d(k,33)  d{k,34)  d(k,35) 

d (k, 36)  d (k, 37 )  d(k,38)  d(k,39)  d(k,40) 

d (k, 41)  d(k, 42)  d(k,43)  d(k,44)  d(k,45)]; 
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function  [A, B, C, D] =abcd_n (a,b,  c,  d,  k) 

%  ' 

%  Name  :  abcd_n.m 
% 

%  Purpose  :  Function  returns  kth  plant  of  full  set  of  plants. 

%  '  Thrust  and  Stall  Margin  Outputs  available  only  the 

%  Fuel  Flow  and  Nozzle  Area  are  available  for  input. 

% 

%  Inputs  :  Must  supply  system  matrices  as  large  matrices  where 
%  each  plant  makes  up  one  row.  User  must  also  supply 

%  which  plant  they  are  interested  in. 

% 

%  Outputs  :  The  desired  plant  in  State  Space  Form 
% 

%  Written  by  :  Zane  D.  Gastineau 


% ================ 


A  =  [a(k,l)  a(k,2)  a(k,3) 
a (k, 4 )  a (k, 5)  a (k, 6) 
a  (k,  7 )  a  (k,  8 )  a(k,9)]; 


B  =  [  b  (k,  1 )  b  (k,  5) 
b (k, 6)  b (k,  10) 
b(k, 11)  b (k, 15)  ] ; 


C  =  [c (k, 1) 
c (k, 4 ) 


c (k, 2)  c (k, 3) 
c (k, 5)  c (k, 6) ]  ; 


D  = 


[d(k, 1)  d(k,  5) 
d (k, 6)  d (k, 10)  ]  ; 
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function  [act, sen] =actblck 


%=======:========:==  = 

% 

%  Name  :  Actblck.m 


%  Purpose  :  To  create  state  space  model  of  the  actuator  block 
%  and  the  sensor  block  for  time  domain  analysis.  All 

%  computational  delays  are  included. 

% 

%  Inputs  :  None  supplied. 

% 

%  Outputs  :  State  Space  models  of  the  actuators  and  sensors. 

% 

%  Written  by  :  Zane  D.  Gastineau 
%  6906  Custer  Rd  #2104 

%  Plano,  TX  75023 

%  ( 937 ) 4  91-377  6 

% 


%  Actuator  Model  #1 
%  Fuel  Flow  Model 

%===========================================================- 

numl= [18] ; 
denl= [1  18]; 

% 

%  Total  60msec  delay  in  this  channel 
%  40  msec  for  bypass  dynamics  and  combustion  delay 

%  20  msec  for  computational  delay 

o. 

o 

[num2, den2] =pade (0 . 02,  4 ) ; %  4th  order  approx,  to  20msec  delay 
[num,  den]  =pade (0.04, 4) ;%  4th  order  approx,  to  4  0msec  delay 
[all,bll, ell, dll]=tf2ss (numl, deni) ; 

[a22, b22, c22, d22] =tf2ss (num, den) ; 

[a33, b33, c33, d33] =tf2ss (num2, den2) ; 

% 

%  Create  state  space  actuator  model 
% 

[al, bl, cl, dl] =series (a33, b33, c33, d33, all , bll, ell, dll) ; 

[al, bl, cl, dl] =se ries (al, bl, cl, dl , a22, b22, c22, d22) ; 


% =================== 

%  Actuator  Model  #2 
%  A8  Actuator 


numl= [ 18 ]  ; 
denl=[l  18]; 

% 

%  Computational  time  delay 
%  20  milliseconds 

% 

[num, den] =pade (0 . 02, 4 ) ;  %4th  order  approx,  to  20msec  delay 
[all,bll,cll,dll]=tf2ss (numl, deni) ; 

[a22,b22, c22, d22] =tf2ss (num, den) ; 
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%  Create  state  space  model  of  actuator 
% 

[a2,b2, c2,d2]=series  (a22, b22, c22, d22, all, bll, ell, dll) ; 

%====  =======:======:=====  =======  =  ====  ============:====== 

%  Sensor  Model 

%  Used  for  both  Thrust  and  Stall  Margin 
%  Feedback.  A  time  delay  is  used  to  represent 
%  the  computational  time  delay  of  the  computer 
%  to  compute  Thrust  and  Stall  Margin.  This 
%  value  is  set  to  20  milliseconds . 


num= [25] ; 
den=[l  25]; 

[num2, den2] =pade (0.02,4) ;  %4th  order  approx,  to  20  msec  delay 
[al,bl, cl, dl] =tf2ss (num,den) ; 

[as,bs, cs,ds]=tf2ss (num2,den2)  ; 

[as,bs, cs, ds] =series  (al, bl, cl, dl, as, bs, cs, ds ) ; 

% 

%  Now  make  state  space  model  for  the  actuator /delay  block: 

% 

[Aact, Bact, Cact, Dact ] =append (al, bl, cl, dl, a2,b2, c2, d2) ; 

% 

%  Now  make  state  space  model  for  the  sensors 
% 

[Asen,  Bsen,  Csen, Dsen] -append (as, bs, cs, ds, as,bs, cs,  ds) ; 


%  Return  Matrices  in  Compact  Form.  Routines  in  Robust  Control 
%  Toolbox . 

act=mksys (Aact , Bact , Cact , Dact ) ; 
sen=mksys (Asen, Bsen, Csen, Dsen) ; 


function  [act, sen] =actblck_w (w) 


%  Name  :  Actblck_w.m 
% 

%  Purpose  :  To  return  complex  evaluation  of  acturators  and 
%  sensors  for  frequency  domain  analysis.  All 

%  computational  delays  are  included. 


%  Inputs  :  User  supplied  frequency  value  in  radians/second. 

%  The  frequency  value  must  be  a  scalar. 

% 

%  Outputs  :  Complex  number  evaluations  for  the  actuators  and 
%  sensors  in  matrix  form. 

% 

%  Written  by  :  Zane  D.  Gastineau 
%  6906  Custer  Rd  #2104 

%  Plano,  TX  75023 

%  (937)491-3776 


% 

%========- 


%  Hous keeping 


s-sqrt (-1) *abs (w) ;  %  if  w  is  complex 

[nr, nc] =size (w) ; 
if  nr  >  1 | nc  >  1 

disp ( 1  frequency  value  must  be  a  scalar’) 
else 

%  Actuator  Model  #1 
%  Fuel  Flow  Model 
% 

%  Total  60msec  delay  in  this  channel 
%  40  msec  for  bypass  dynamics  and  combustion  delay 

%  20  msec  for  computational  delay 

actl= ( 18 . / (s  +  18 ) ) . *exp (-s*60e-3)  ; 

%==========--=-===  ===:=-==-====  =  -  =  -===========:========: 

%  Actuator  Model  #2 
%  A8  Actuator 
% 

%  Computational  time  delay 
%  20  milliseconds 

act2=(18./ (s+18) ) . *exp (-s*20e-3) ; 
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%  Sensor  Model 

%  Used  for  both  Thrust  and  Stall  Margin 
%  Feedback.  A  time  delay  is  used  to  represent 
%  the  computational  time  delay  of  the  computer 
%  to  compute  Thrust  and  Stall  Margin.  This 
%  value  is  set  to  20  milliseconds. 


senl=exp (-s*20e-3) ; 

act=diag ( [actl  act2] ) ; 
sen=diag ( [senl  senl] ) ; 
end 


121 


function  [cont, decple] =gk 


% 

%  Name  :  gk.m 
% 

%  Purpose  :  Routine  to  return  state  space  form  of  controller 
%  and  decoupler 

% 

%  Inputs  :  No  inputs  required 
.% 

%  Outputs  :  State  space  models  of  the  precompensator  and 
%  controller. 

% 

%  Written  by  :  Zane  D.  Gastineau 


%  Diagonal  Controller 

gllnum=7*conv ( [1/100  0.15  1] , conv ( [1/4 . 201  1] , [1/4.958  1] ) ) ; 

glldena=conv ( [1  0],conv([l/5  1] , conv {[ 1/8 . 560  1], [1/69.06  1 ] ) ) ) ; 

gllden=conv (glldena, [1/120  1] ) ; 

g22num=0 . 8*conv ([1/0.7025  1] , [1/2.38  1]); 

g22den-conv( [1  0] , conv ( [1/0 . 5124  1], [1/6.723  1] ) ) ; 

[glia, glib, gllc, gild] =tf2ss (gllnum, gllden) ; 

[g22a, g22b, g22c, g22d] -tf2ss (g22num, g22den) ; 

[Ag, Bg, Cg, Dg] ^append (glia, glib, gllc, gild, g22a,  g22b, g22c,  g22d) ; 

%  Precompensator 

kllnum=l; 
kllden=l ; 
k22num=l ; 
k22den=l; 

kl2num=0 . 05*conv ( [1/0.1034  1],  [1/2.144  1]  )  ; 
kl2den=conv( [1/1.547  1], [1/29.86  1 ] ) ; 
k21num=0. 005* [1/0.02  1] ; 
k21den= [1/2 . 087  1] ; 

[klla, kllb, kllc, klld] =tf 2ss ( kllnum, kllden) ; 

[k22a, k22b, k22c, k22d] =tf2ss (k22num, k22den) ; 

[kl2a, kl2b, kl2c, kl2d] =tf2ss (kl2num, kl2den) ; 

[k21a, k21b, k21c, k21d] =tf2ss (k21num, k21den) ; 

skl=size (klla) ; sk2=size (kl2a) ; sk3=size (k21a) ; sk4=size (k22a) ; 
Ako=[kl2a  zeros (sk2 ( 1) , sk3 (2 ) ) ; 

zeros (sk3 (1) , sk2 (2) )  k21a] ; 

Bko= [zeros (sk2 (1) , 1)  kl2b 
k21b  zeros (sk3 (1)  ,  1)  ]  ; 

Cko=[kl2c  zeros (1, sk3 (2) ) 
zeros (1, sk2 (2) )  k21c] ; 

Dko= [ 0  kl2d; k21d  0] ; 

[Akd, Bkd, Ckd, Dkd] =append (klla, kllb, kllc, klld, k22a, k22b, k22c, k22d) ; 
[Ak, Bk, Ck, Dk] =parallel (Akd, Bkd, Ckd, Dkd,  Ako,  Bko,  Cko,  Dko)  ; 
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%  Return  in  compact  form 
%  Routines  in  Robust  control  toolbox 


cont=mksys (Ag, Bg, Cg,  Dg) ; 
decple=mksys (Ak,  Bk,  Ck,  Dk) ; 
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function  [Ctrl, dec] =gk_w (w) 


% 

%  Name  :  gk_w.m 
% 

%  Purpose  :  Routine  to  return  frequency  response  of  controller 
%  and  decoupler 

% 

%  Inputs  :  supply  desired  frequency  to  compute  at 
% 

%  Outputs  :  Complex  matrix  representing  the  precompensator  and 
%  Controller. 

% 

%  Written  by  :  Zane  D.  Gastineau 
% 

%========================:==================::======:=:===============:==========:=:: 

%  Housekeeping 
s=sqrt (-1) . *w; 

%================================================= 

%  Controller  responses 

%===^========-=“===“==================“=“=============:==“=! 


gllnuma=7 . * ( (s/100) .  A2+ (2*0 . 75/10) . *s+l) ; 
gllnumb- (s/4 .201+1) .* (s/4 . 958+1) ; 
glldena=s . * (s/5+1) .*(s/8.560+l)  ; 
glldenb- (s/69. 06+1) .* (s/120+1) ; 

gll=(gllnuma. *gllnumb) . / (glldena . *glldenb)  ; 

g22num=0.8.* (s/0.7025+1) .* (s/2 . 38+1 ) ; 
g22den-s.* (s/0.5124+1) .* (s/6 . 723+1 )  ; 

g22=g22num. /g22den; 

%  Precompensator  Responses 


kll=l; k22=l; 


k!2num=0 .05.*  (s/0.1034  +  1)  .*  (s/2 . 14  4  +  1)  ;■ 
kl2den= (s/1.547  +  1) . * (s/2  9 . 86+1 )  ; 

kl2=kl2num. /kl2den; 

k21num=0. 005.*  (s/0.02+1) ; 
k21den= (s/2.087+1)  ; 

k21=k21num. /k21den; 
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%  returned  matrices 


ctrl= [gll  0; 0  g22] ; 
dec= [kll  k!2 ; k21  k22] ; 
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function  [F]=pfilter 
% 

%  Name  :  pfilter.m 
% 

%  Purpose  :  Routine  to  return  state  space  model  of  the 
%  prefilter 

% 

%  Inputs  :  No  inputs  required. 

o 

"o 

%  Outputs  :  State  space  models  of  the  prefilter. 

% 

%  Written  by  :  Zane  D.  Gastineau 
f llnum=l; 

fllden= [1/3 . 963  1]  ; 

% 

f22num=conv( [1/2.703  1], [1/11.96  1]); 
f22den=conv( [1/1. 32  1], [1/3.139  1]); 

% 

[f 11a, f lib, flic, f lid] ~tf2ss {f llnum, fllden) ; 

[f 22a, f22b, f22c, f22d] =tf2ss (f22num, f22den) ; 

% 

[Fa,  Fb,  Fc, Fd] =append ( f 11a, f lib, f 11c, f lid, f 22a, f22b, f22c, f22d) ; 


%=====“==““==“=====s===“===“==“: 
%  Return  in  Compact  form 
%  Routines  in  Robust  control  toolbox 

F=mksys (Fa, Fb, Fc, Fd) ; 
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function  [F] =pf ilter_w (w) 

% 

%  Name  :  pfilter_w 
% 

%  Purpose  :  Routine  to  return  frequency  response  of  the 
%  prefilter. 

% 

%  Inputs  :  Supply  desired  frequency  to  compute  at. 

% 

%  Outputs  :  Complex  matrix  representing  the  prefilter. 

% 

%  Written  by  :  Zane  D.  Gastineau 

%=-========!=======-=====:=====================:========:==============:== 

%  Housekeeping 

%  =~~~~=====— ~~=====——==—— —————— ====—~——=~==—~~~~~~~~~~::=:~== 

s=sqrt (-1) . *w; 

%=======^======  ===========-==:=:=====-==-===  ====---=========:====-==========: 

%  Prefilter  responses 

%-======================:======================================== 

f llnum=l; 

fllden=(s/3. 963+1) ; 

f 11-fllnum. /f llden; 

f22num= (s/2. 703+1) . * (s/11.96+1)  ; 
f22den=( s/1. 32  +  1) .*  (s/3. 139+1) ; 

f22=f22num. /f22den; 

%  Returned  Matrices 
F= [f 11  0 ; 0  f22] ; 
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<A°  C*° 


function  [lambda,  v]=perron(Z) 


o 


Name  :  Perron. m 

%  Purpose  :  Extracts  the  Perron  root  if  it  exists 
% 

%  Inputs  :  Complex  Matrix 
% 

%  Outputs  :  The  Perron  root 
% 

Z=abs (Z) ; 

if  (nargout==l) 

lambda=real (max (eig (Z) ) ) ; 

else 

[X, L] =eig (Z) ; 

L=diag (L) ; 

[L,  index] =sort (L) ; 
lambda=re£l (L (length (L) ) ) ; 
x=real (X( : , index (length (L) ))); 

[Y, L] =eig ( Z 1 ) ; 

L=diag (L) ; 

[L, index] =sort (L) ; 

y=real (Y ( : , index (length (L) ) ) ) ; 

v= (x*y ’ ) / (y 1 *x) ; 

end 
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APPENDIX  D 


NORMALIZED  PLANT  MODELS 

The  following  pages  contain  the  normalized  linear  model  data  for  the  turbojet 
engine  used  for  this  work.  This  data  was  produced  from  a  nonlinear  simulation  built  in 
the  Math  Works  Simulink  simulation  environment  using  the  linearization  tools  that  were 
available.  All  the  points  represent  the  engine  at  sea  level  static  conditions  and  the  points 
are  generated  at  values  of  corrected  speed  from  idle  (70%  corrected  speed)  to  full  speed 
(100%  corrected  speed).  These  models  do  not  include  the  actuators  and  sensors.  The 
states,  inputs  and  outputs  are  given  below. 

Rotational  Speed  (N) 

Combustor  Pressure  (P4 ) 

Tailpipe  Pressure  (P6) 

Fuel  Flow  (wf) 

Nozzle  Area  (Ag) 

Thrust  (Fn) 

Stall  Margin  (SM) 

Rotational  Speed  (N) 

Compressor  Exit  Temp.  (T3) 

Turbine  Exit  Temp.  (T5) 

Turbine  Exit  Press.  (P5) 

Compressor  Exit  Press.  (P? ) 

Turbine  Inlet  Temp.  (T4) 

%  Corrected  Speed  (%Ncorr) 
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0.0098  -0.0051  1.6550  0.0119  0.0426 

2.1684  -2.8139  0.0000  0.0000  0.0000 

1.0000  0.0000  0.0000  0.0000  0.0000 
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-0.4802  0.7523  0.0000  4.2205  0.0000 

0.6044  0.0000  0.0000  0.0000  0.0000 


0.0121  -0.0097  1.3644  0.0154  0.0781 

1.8082  -3.0284  0.0000  0.0000  0.0000 

1.0000  0.0000  0.0000  0.0000  0.0000 

•3.0797  6.7650  -10.7836  2.0685  0.0000  0.1594  0.4314  0.0000  0.0000  0.0000 
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0.6044  0.0000  0.0000  0.0000  0.0000 


0.0164  -0.0124  1.2323  0.0179  0.1142 

1.8587  -2.8451  0.0000  0.0000  0.0000 

1.0000  0.0000  0.0000  0.0000  0.0000 
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0.6044  0.0000  0.0000  0.0000  0.0000 


0.0149  -0.0120  0.5509  0.0299  0.5282 

0.2109  -0.8924  0.0000  0.0000  0.0000 

1.0000  0.0000  0.0000  0.0000  0.0000 
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NOMENCLATURE 


A  Coefficient  System  Matrix 

A  Area  (ft2) 

Aij(co)  Lower  Bound  of  System  Specification 

B  Input  System  Matrix 

By(co)  Upper  Bound  of  System  Specification 

C  Output  System  Matrix 

C  Controller,  Feedback  Regulator 

C  Proportionality  Constant 

D  Transmission  System  Matrix 

DoD  Department  of  Defense 

F  Prefilter  Matrix 

Fn  Thrust 

G  Diagonal  Controller  Matrix 

G  Proportional  Gain  Matrix  for  Proportional-Plus-Integral  Observer 
IEC  Intelligent  Engine  Control 

IHPTET  Integrated  High  Performance  Turbine  Engine  Technology 

IMC  Internal  Model  Control 

INA  Inverse  Nyquist  Array 

J  Moment  of  Inertia 
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K 

K 

K 

M 

MBC 

MRAC 

N 

N 

NASA 

P 

P 

PM 

A 

P 

P 

Q 

QFT 

QNA 

R 

S 

SAE 

SM 

T 

T 


Precompensating  Matrix 
Proportionality  Constant 

Integral  Gain  Matrix  for  Proportional-Plus-Integral  Observer 

Reference  Model 

Model-Based  Control 

Model-Reference  Adaptive  Control 

Speed  (rpm) 

Acceleration  (rpm/sec) 

National  Aeronautics  and  Space  Administration 
Plant  Matrix 
Pressure  (lb/in2) 

Plant  Model 

Inverse  of  the  Matrix  P 

Family  of  Plants,  Plant  Set 

Heating  Value  of  Fuel 

Quantitative  Feedback  Theory 

Quantitative  Nyquist  Array 

Scaling  Matrix 

Diagonal  Sensitivity  Matrix 

Society  of  Automotive  Engineers 

Stall  Margin 

Temperature 

Diagonal  Complimentary  Sensitivity  Function 
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T  Transfer  Function  Matrix 

V  Velocity 

V  Volume 

cp  Specific  Heat  at  Constant  Pressure 
e  Error 
h  Enthalpy 

k  Ratio  of  Specific  Heats 
i{  Loop  Transmission  Element  for  the  ith  Loop 
m  Mass  flow  (lbm/sec) 
qy  ith  row,  jth  column  of  P  "1 
s  Laplace  Operator 

u  Input  Vector 

w  Flow 

x  State  Variable  Vector 

x  First  Derivative  of  State  Variable  V  ector 
x  Estimated  State  Variables 
y  Output  Vector 

a  Vector  of  Parameters  Corresponding  to  Points  of  Linearization 
y  Interaction  Index 

r)  Efficiency 

X  Eigenvalue 

A,  Maximum  Eigenvalue 


143 


A,p  Perron  Root 

p  Spectral  Radius 
t  Torque 

0  frequency  (rad/sec) 

Subscripts 

b  Burner,  Combustor 
c  Compressor 
d  Diagonal 
f  Fuel 
i  ith  element 

ij  ith  input,  jth  output 

o  Off  Diagonal 

t  Turbine 

Y/R  Output  to  Reference 
0-9  Engine  Station  Designations 
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